Making Generic Commands

Let’s let put the Bash shell to work for us by being smart with our variables

Shell Variables

Assign the variables in this notebook.

source bioinf_intro_config.sh
mkdir -p $TRIMMED $STAR_OUT

Using a FASTQ variable

Most of the commands in our pipeline are complicated! To apply our pipeline to different FASTQs, we need to change things in multiple places. For example, just to run trimming with fastq-mcf, we need to change things in two places between each run: the input FASTQ and the output FASTQ. If we were doing this with paired-end reads, we would have to change four things. Doing this by hand is not only tedious, but error prone. Doing almost the same thing repeatedly is something that people are bad at, but computers are very good at! So let’s get the computer to do the hard work. Because the Bash shell is almost magical (it is a full fledged programming language), we can do this easily.

FASTQ="A"
echo "RUNNING FASTQ: ~~~~${FASTQ}~~~~"
echo "TRIMING: ${FASTQ}"
echo "MAPPING: ${FASTQ}"

RUNNING FASTQ: ~~~~A~~~~
TRIMING: A
MAPPING: A

The for loop in Bash is conceptually the same as in any other programming language, although the syntax may be different. The do and done are essential - do needs to be before the “loop body” (what is going to be repeated) and done needs to be after it.

So let’s try something almost useful:

FASTQ="21_2019_P_M1_S21_L002_R1"
echo "RUNNING FASTQ: ~~~~${FASTQ}~~~~"
echo "TRIMING: ${FASTQ}"
echo "MAPPING: ${FASTQ}"

RUNNING FASTQ: ~~~~21_2019_P_M1_S21_L002_R1~~~~
TRIMING: 21_2019_P_M1_S21_L002_R1
MAPPING: 21_2019_P_M1_S21_L002_R1

Now for the real thing …

Let’s run fastq-mcf:

FASTQ="21_2019_P_M1_S21_L002_R1"
echo "TRIMING: $FASTQ"
fastq-mcf $MYINFO/neb_e7600_adapters.fasta \
    $RAW_FASTQS/${FASTQ}_001.fastq.gz \
        -q 20 -x 0.5 \
        -o $TRIMMED/${FASTQ}_001.trim.fastq.gz

TRIMING: 21_2019_P_M1_S21_L002_R1
Command Line: /home/jovyan/work/scratch/bioinf_intro/myinfo/neb_e7600_adapters.fasta /data/hts_2019_data/hts2019_pilot_rawdata/21_2019_P_M1_S21_L002_R1_001.fastq.gz -q 20 -x 0.5 -o /home/jovyan/work/scratch/bioinf_intro/trimmed_fastqs/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz
Scale used: 2.2
Phred: 33
Threshold used: 751 out of 300000
Adapter Adapter (AGATCGGAAGAGCACACGTCTGAACTCCAGTCA): counted 2515 at the 'end' of '/data/hts_2019_data/hts2019_pilot_rawdata/21_2019_P_M1_S21_L002_R1_001.fastq.gz', clip set to 6
Files: 1
Total reads: 2437108
Too short after clip: 1347
Clipped 'end' reads: Count: 44977, Mean: 15.55, Sd: 8.27
Trimmed 288960 reads by an average of 1.70 bases on quality < 20

Now let’s do the same thing for STAR

FASTQ="21_2019_P_M1_S21_L002_R1"
echo "MAPPING: $FASTQ"
STAR \
    --runMode alignReads \
    --twopassMode None \
    --genomeDir $GENOME_DIR \
    --readFilesIn $TRIMMED/${FASTQ}_001.trim.fastq.gz \
    --readFilesCommand gunzip -c \
    --outFileNamePrefix ${STAR_OUT}/${FASTQ}_ \
    --quantMode GeneCounts \
    --outSAMtype None \
    --runThreadN 2

MAPPING: 21_2019_P_M1_S21_L002_R1
Jun 26 16:18:39 ..... started STAR run
Jun 26 16:18:39 ..... loading genome
Jun 26 16:18:39 ..... started mapping
Jun 26 16:20:03 ..... finished successfully

And let’s check the result

ls ${STAR_OUT}

21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L001_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L001_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L001_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_Aligned.out.bam
21_2019_P_M1_S21_L002_R1_Log.final.out
21_2019_P_M1_S21_L002_R1_Log.out
21_2019_P_M1_S21_L002_R1_Log.progress.out
21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_SJ.out.tab
genome_Log.out
multiqc_data
multiqc_report.html

head ${STAR_OUT}/${FASTQ}_ReadsPerGene.out.tab

N_unmapped      46060   46060   46060
N_multimapping  35006   35006   35006
N_noFeature     12466   2145291 18783
N_ambiguous     203327  820     316
CNAG_04548      0       0       0
CNAG_07303      0       0       0
CNAG_07304      6       0       6
CNAG_00001      0       0       0
CNAG_07305      1       0       1
CNAG_00002      51      0       51