== Creating an analysis pipeline of compressed files ==

Almost all steps of a short-read analysis pipeline can be created with gzipped files to save storage space.  For example,

# Untar from solexa_public, gunzip, and gzip \\
tar -xOzf /lab/solexa_public/LAB/RUN/QualityScore/Foo.1.tar.gz | gzip -f > fastq/My_sample_2.1.fq.gz \\
tar -xOzf /lab/solexa_public/LAB/RUN/QualityScore/Foo.2.tar.gz | gzip -f > fastq/My_sample_2.2.fq.gz \\

# Run FastQC \\
fastqc fastq/My_sample_2.1.fq.gz \\
fastqc fastq/My_sample_2.2.fq.gz \\

# Trim adapters \\
gunzip -c fastq/My_sample_2.1.fq.gz | fastx_clipper -v -z -l 30 -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG -o fastq/My_sample_2.noVec.1.fq.gz \\
gunzip -c fastq/My_sample_2.2.fq.gz | fastx_clipper -v -z -l 30 -a GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG -o fastq/My_sample_2.noVec.2.fq.gz \\
fastqc fastq/My_sample_2.noVec.1.fq.gz \\
fastqc fastq/My_sample_2.noVec.2.fq.gz \\

# Drop reads of low quality \\
gunzip -c fastq/My_sample_2.noVec.1.fq.gz | fastq_quality_filter -v -q 20 -p 75 -z -o fastq/My_sample_2.noVec.qgt20.1.fq.gz \\
gunzip -c fastq/My_sample_2.noVec.2.fq.gz | fastq_quality_filter -v -q 20 -p 75 -z -o fastq/My_sample_2.noVec.qgt20.2.fq.gz \\

# Trim low-quality bases at end of reads \\
gunzip -c fastq/My_sample_2.noVec.qgt20.1.fq.gz | fastq_quality_trimmer -v -t 20 -l 30 -z -o fastq/My_sample_2.noVec.qgt20.t.1.fq.gz \\
gunzip -c fastq/My_sample_2.noVec.qgt20.2.fq.gz | fastq_quality_trimmer -v -t 20 -l 30 -z -o fastq/My_sample_2.noVec.qgt20.t.2.fq.gz \\

# After trimming, get reads that are still paired (and output as regular fastq) \\
/nfs/BaRC_Public/BaRC_code/Perl/cmpfastq/cmpfastqgz.pl fastq/My_sample_2.noVec.qgt20.t.1.fq.gz fastq/My_sample_2.noVec.qgt20.t.2.fq.gz

# Map paired reads with bowtie \\
# THIS STEP CANNOT USE GZ FILES \\
bowtie -l 30 -n 1 -X 1000 --best -S /nfs/genomes/a.thaliana_TAIR_10/bowtie/tair10_SI -1 fastq/My_sample_2.noVec.qgt20.t.1.fq-common.out -2 fastq/My_sample_2.noVec.qgt20.t.2.fq-common.out mapped_reads/My_sample_2.30.1.sam

# Sort and index mapped reads \\
/nfs/BaRC_Public/BaRC_code/Perl/SAM_to_BAM_sort_index/SAM_to_BAM_sort_index.pl mapped_reads/My_sample_2.30.1.sam \\
# Delete SAM file (since BAM file has all the same info) \\
rm -f mapped_reads/My_sample_2.30.1.sam \\

# Get summary counts \\
zcat fastq/My_sample_2.1.fq.gz | echo $((`wc -l`/4)) \\
zcat fastq/My_sample_2.noVec.1.fq.gz | echo $((`wc -l`/4)) \\
zcat fastq/My_sample_2.noVec.2.fq.gz | echo $((`wc -l`/4)) \\
zcat fastq/My_sample_2.noVec.qgt20.1.fq.gz | echo $((`wc -l`/4)) \\
zcat fastq/My_sample_2.noVec.qgt20.2.fq.gz | echo $((`wc -l`/4)) \\
more fastq/My_sample_2.noVec.qgt20.t.1.fq | echo $((`wc -l`/4)) \\
more fastq/My_sample_2.noVec.qgt20.t.2.fq | echo $((`wc -l`/4)) \\
zmore fastq/My_sample_2.noVec.qgt20.t.1.fq-common.out.gz | echo $((`wc -l`/4)) \\
samtools flagstat mapped_reads/My_sample_2.30.1.sorted.bam > mapped_reads/My_sample_2.30.1.flagstat.txt \\
grep ' mapped (' mapped_reads/My_sample_2.30.1.flagstat.txt