Matt:LabNotes/2015-4-8
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
- Output:
CountReadsPer_Gene_Probe_excludeSoftMask.pl
- Output:
- CA12kNov2014_V7_SoftMaskedProbes.txt
- All probes containing softmasked regions, '1' in 2nd column indicates it also is a part of CA12kNov2014_V4_SoftMaskedProbes.txt
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_Excluded_Probecounts.txt
- Counts for only probes that are not in CA12kNov2014_V7_SoftMaskedProbes.txt
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_Excluded_Genecounts.txt
- Counts for only genes that have no SoftMaskedProbes
- CA12kNov2014_V7gDNA_R1_H1H2_sorted_filtered_ExcludedProbesOnly_Genecounts
- Counts for all genes but do not include count for SoftMaskedProbes
- CA12kNov2014_V7_SoftMaskedProbes.txt