Matt:LabNotes/2015-4-8

From ZhangLabWiki
Revision as of 20:20, 8 April 2015 by >Mzcai (Created page with "==CA12k_Nov2014_V7 in vitro Capture Sequencing Analysis== *With in vitro Capture data, DARTFISH gene counts can be normalized and compared to other methods such as RNA-Seq *Fo...")
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

CA12k_Nov2014_V7 in vitro Capture Sequencing Analysis

  • With in vitro Capture data, DARTFISH gene counts can be normalized and compared to other methods such as RNA-Seq
  • For V4 in vitro Capture Sequencing Analysis see: Matt:LabNotes/2015-3-19

Mapping MiSeq reads to Oligo Sequences

Convert Probelist to Fasta File

CA12kNov2014_Probelist2Fasta.pl

 bowtie2-build CA12k_Nov2014_V7_H1H2.fa CA12k_Nov2014_V7_H1H2

Map with Bowtie2

Mapping full read to reference
  • Had a 0% mapping rate
    • Reason is because the tail end of 100bp reads had very low base quality scores (some reads had up 80 "#")
    • Forums call this "B-tailing" because in Phred64 'B' represents score of 2 while Phred33 uses '#'
    • This means something went wrong with sequencing at a certain cycle (possibly due to fragment being shorter than read length but this library should definitely be long enough)
Use fastx toolkit to visualize read quality stats
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_quality_stats/fastx_quality_stats -Q33 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA_qualstats.txt
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_quality_stats/fastx_quality_stats -Q33 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA_qualstats.txt
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/scripts/fastq_quality_boxplot_graph.sh -i MC-20150121_CA12kNov2014_V7gDNA_qualstats.txt -o MC-20150121_CA12kNov2014_V7gDNA_qualstats.png -t CA12kNov2014_V7_gDNA
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/scripts/fastq_quality_boxplot_graph.sh -i MC-20150121_CA12kNov2014_V7cDNA_qualstats.txt -o MC-20150121_CA12kNov2014_V7cDNA_qualstats.png -t CA12kNov2014_V7_cDNA

File:MC-20150121 CA12kNov2014 V7gDNA qualstats.png File:MC-20150121 CA12kNov2014 V7cDNA qualstats.png

Trim reads with fastx toolkit
  • trimmed.fastq: Trim by base quality (minimum quality=10, minimum length=20)
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastq_quality_trimmer/fastq_quality_trimmer -Q33 -v -t 10 -l 20 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmed.fastq
 Minimum Quality Threshold: 10
 Minimum Length: 20
 Input: 2390629 reads.
 Output: 2390528 reads.
 discarded 101 (0%) too-short reads.
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastq_quality_trimmer/fastq_quality_trimmer -Q33 -v -t 10 -l 20 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq
 Minimum Quality Threshold: 10
 Minimum Length: 20
 Input: 3675652 reads.
 Output: 3675504 reads.
 discarded 148 (0%) too-short reads.
  • trimmedv2.fastq: Trim by fixed length (keep bases 2-45; length = 44)
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_trimmer/fastx_trimmer -Q33 -f 2 -l 45 -i MC-20150121_CA12kNov2014_V7gDNA-4.R1.fastq -o MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq
 /home/kunzhang/softwares/fastx_toolkit-0.0.13.2/src/fastx_trimmer/fastx_trimmer -Q33 -f 2 -l 45 -i MC-20150121_CA12kNov2014_V7cDNA-5.R1.fastq -o MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq
Map quality trimmed reads to reference
  • Mapping rate still too low
    • Looking at fastq by eye shows that the algorithm is missing easy ones that I can find by ctrl-f
 bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmed.fastq > CA12kNov2014_V7gDNA_R1_H1H2.sam 2> CA12kNov2014_V7gDNA_stderr.txt &
 2390528 reads; of these:
 2390528 (100.00%) were unpaired; of these:
   2199579 (92.01%) aligned 0 times
   190821 (7.98%) aligned exactly 1 time
   128 (0.01%) aligned >1 times
 7.99% overall alignment rate
 bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq > CA12kNov2014_V7cDNA_R1_H1H2.sam 2> CA12kNov2014_V7cDNA_stderr.txt &
 3675504 reads; of these:
 3675504 (100.00%) were unpaired; of these:
   3366177 (91.58%) aligned 0 times
   309126 (8.41%) aligned exactly 1 time
   201 (0.01%) aligned >1 times
 8.42% overall alignment rate
Map v2 trimmed reads to reference
  • Much better mapping rate but examining fastq by eye still shows the algorithm can do better
    • A common one missed is PDE1A because it has a repeating "AC" (softmasked) that I guess has a lot of sequencing error
 bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq > CA12kNov2014_V7gDNA_R1_H1H2.sam 2> CA12kNov2014_V7gDNA_stderr.txt &
 2390629 reads; of these:
 2390629 (100.00%) were unpaired; of these:
   191725 (8.02%) aligned 0 times
   2189100 (91.57%) aligned exactly 1 time
   9804 (0.41%) aligned >1 times
 91.98% overall alignment rate
 bowtie2 --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq > CA12kNov2014_V7cDNA_R1_H1H2.sam 2> CA12kNov2014_V7cDNA_stderr.txt &
 3675652 reads; of these:
 3675652 (100.00%) were unpaired; of these:
   275464 (7.49%) aligned 0 times
   3385203 (92.10%) aligned exactly 1 time
   14985 (0.41%) aligned >1 times
 92.51% overall alignment rate
Map v2 trimmed reads with Local alignment
  • Best mapping rate so far, I guess local alignment helps with all the sequencing errors in the reads?
    • Overall alignment rate comparable to V4 in vitro capture alignment rate
    • High multiple alignments but algorithm should choose the best mapping score
  • USING THIS SAM FILE FOR FURTHER ANALYSIS
 bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7gDNA-4.R1.trimmedv2.fastq > CA12kNov2014_V7gDNA_R1_H1H2_localv2.sam 2> CA12kNov2014_V7gDNA_stderr_localv2.txt &
 2390629 reads; of these:
 2390629 (100.00%) were unpaired; of these:
   76473 (3.20%) aligned 0 times
   886116 (37.07%) aligned exactly 1 time
   1428040 (59.73%) aligned >1 times
 96.80% overall alignment rate
 bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmedv2.fastq > CA12kNov2014_V7cDNA_R1_H1H2_localv2.sam 2> CA12kNov2014_V7cDNA_stderr_localv2.txt &
 3675652 reads; of these:
 3675652 (100.00%) were unpaired; of these:
   109047 (2.97%) aligned 0 times
   1368908 (37.24%) aligned exactly 1 time
   2197697 (59.79%) aligned >1 times
 97.03% overall alignment rate
 samtools view -bS CA12kNov2014_V7gDNA_R1_H1H2_localv2.sam | samtools sort - CA12kNov2014_V7gDNA_R1_H1H2_sorted
 samtools view -h -F 4 CA12kNov2014_V7gDNA_R1_H1H2_sorted.bam > CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered.sam
 samtools view -bS CA12kNov2014_V7cDNA_R1_H1H2_localv2.sam | samtools sort - CA12kNov2014_V7cDNA_R1_H1H2_sorted
 samtools view -h -F 4 CA12kNov2014_V7cDNA_R1_H1H2_sorted.bam > CA12kNov2014_V7cDNA_R1_H1H2_sorted_filtered.sam
  • Using local alignment for quality trimmed reads leads to too many multiple alignments
 bowtie2 --local --phred33 -x CA12k_Nov2014_V7_H1H2 -q MC-20150121_CA12kNov2014_V7cDNA-5.R1.trimmed.fastq > CA12kNov2014_V7cDNA_R1_H1H2_local.sam 2> CA12kNov2014_V7cDNA_stderr_local.txt &
 3675504 reads; of these:
 3675504 (100.00%) were unpaired; of these:
   11255 (0.31%) aligned 0 times
   64099 (1.74%) aligned exactly 1 time
   3600150 (97.95%) aligned >1 times
 99.69% overall alignment rate

Count # of Reads Aligned to each Probe

CountReadsPer_Gene_Probe.pl

CountReadsPer_Gene_Probe_excludeSoftMask.pl

Regression Analsis