Working with Loops

Let’s kick it up another notch we have lots of FASTQs, let’s run our analysis on more than one!

Shell Variables

Assign the variables in this notebook.

source bioinf_intro_config.sh
mkdir -p $TRIMMED $STAR_OUT

for FASTQ in 21_2019_P_M1_S21_L002_R1 21_2019_P_M1_S21_L003_R1
    do
        echo "RUNNING FASTQ: ${FASTQ}"
    done

RUNNING FASTQ: 21_2019_P_M1_S21_L002_R1
RUNNING FASTQ: 21_2019_P_M1_S21_L003_R1

Now let’s run both samples throught the pipeline:

for FASTQ in 21_2019_P_M1_S21_L002_R1 21_2019_P_M1_S21_L003_R1
    do
        echo "---------------- TRIMMING: $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

        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
    done

---------------- TRIMMING: 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
---------------- MAPPING: 21_2019_P_M1_S21_L002_R1 ----------------
Jun 26 16:26:49 ..... started STAR run
Jun 26 16:26:49 ..... loading genome
Jun 26 16:26:50 ..... started mapping
Jun 26 16:28:13 ..... finished successfully
---------------- TRIMMING: 21_2019_P_M1_S21_L003_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_L003_R1_001.fastq.gz -q 20 -x 0.5 -o /home/jovyan/work/scratch/bioinf_intro/trimmed_fastqs/21_2019_P_M1_S21_L003_R1_001.trim.fastq.gz
Scale used: 2.2
Phred: 33
Threshold used: 751 out of 300000
Adapter Adapter (AGATCGGAAGAGCACACGTCTGAACTCCAGTCA): counted 2552 at the 'end' of '/data/hts_2019_data/hts2019_pilot_rawdata/21_2019_P_M1_S21_L003_R1_001.fastq.gz', clip set to 6
Files: 1
Total reads: 2507049
Too short after clip: 1463
Clipped 'end' reads: Count: 46508, Mean: 15.55, Sd: 8.23
Trimmed 300377 reads by an average of 1.69 bases on quality < 20
---------------- MAPPING: 21_2019_P_M1_S21_L003_R1 ----------------
Jun 26 16:28:29 ..... started STAR run
Jun 26 16:28:29 ..... loading genome
Jun 26 16:28:29 ..... started mapping
Jun 26 16:29:51 ..... 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
21_2019_P_M1_S21_L003_R1_Log.final.out
21_2019_P_M1_S21_L003_R1_Log.out
21_2019_P_M1_S21_L003_R1_Log.progress.out
21_2019_P_M1_S21_L003_R1_ReadsPerGene.out.tab
21_2019_P_M1_S21_L003_R1_SJ.out.tab
genome_Log.out
multiqc_data
multiqc_report.html