Daniel:Notebook/GenomeMiner/2013-9-12: Difference between revisions
>Djacobse |
>Djacobse |
||
(14 intermediate revisions by the same user not shown) | |||
Line 25: | Line 25: | ||
So it looks like even with indisputable base quality and sequences, about 0.2% of sequences are so close that bowtie2 thinks they can be several sequences. | So it looks like even with indisputable base quality and sequences, about 0.2% of sequences are so close that bowtie2 thinks they can be several sequences. | ||
==Bowtie on Perfect Mock Sequence Data([[Daniel:Notebook/GenomeMiner/2013-9-11|Started 09/11/2013]])== | '''Convert Sam to Bam and Sort''' | ||
samtools view -bS mockseq_v4s1_perfseqbqual.sam | samtools sort - mockseq_v4s1_perfseqbqual_sort | |||
'''Filter and Convert Back to Sam''' | |||
samtools view -h -F 4 -q 8 mockseq_v4s1_perfseqbqual_sort.bam > mockseq_v4s1_perfseqbqual_sortfilt.sam | |||
'''Replace Matched Pairs with Equals''' | |||
samtools calmd -eS mockseq_v4s1_perfseqbqual_sortfilt.sam V4S1.fa > mockseqv4s1perfseqbqual_sortfilteq.sam | |||
'''Perl Script Results''' script: [[File:imp_count_mismatch.txt]] | |||
Error Rate: 0 | |||
Error Rate of Insertions: 0 | |||
Error Rate of Deletions: 0 | |||
Error Rate of Substitutions: 0 | |||
Count: 2955000 perfect, 0 errors | |||
==Bowtie on Perfect Mock Sequence Data ([[Daniel:Notebook/GenomeMiner/2013-9-11|Started 09/11/2013]])== | |||
Using Matt's Perlscript: [[File:imp_count_mismatch.txt]] | Using Matt's Perlscript: [[File:imp_count_mismatch.txt]] | ||
Line 33: | Line 50: | ||
Error Rate of Deletions: 0 | Error Rate of Deletions: 0 | ||
Error Rate of Substitutions: 0 | Error Rate of Substitutions: 0 | ||
Counts: 2955000 perfect, 0 with error | |||
So this is expected from the sequencing data with no errors. Interestingly, 4000 probes were not counted at all, meaning (most likely) that 4000 probes were of quality too low (they were filtered out in the calmd stage). 4000 is a suspicious number, being even, which could mean that 4 different probes were not counted. It could also just be coincidence. However, the error is present in reads of good quality and reads of normal quality. | |||
==Bash Script== | |||
To save myself a lot of time as I run through all of this stuff, I've created a bash script to run through the pipeline (above), going through bowtie2 all the way to the sam file with the equals signs. It's here: [[Media:Hl155bash.sh]] | |||
==Analysis on Mimic Data== | |||
Mimic data is named such because it has the same error rates calculated from the actual sequencing data (1.01% substitutions, 0.26% insertions, and 0.05% deletions). It was generated using ''MockHL155_Master, Switch 2'' with parameters 1.01, 0.26, and 0.05. | |||
The analysis was done using the bash script from above. | |||
===Perl Script Results=== | |||
Used Matt's perl script (above) | |||
Error Rate: 0.024866487399162 | |||
Error Rate of Insertions: 3.25182915389907e-05 | |||
Error Rate of Deletions: 1.25070352073041e-06 | |||
Error Rate of Substitutions: 0.0248327184041023 | |||
'''Errors:''' | |||
{| class="wikitable" <hiddentext>generated with [[:de:Wikipedia:Helferlein/VBA-Macro for EXCEL tableconversion]] V1.8</hiddentext> | |||
|- style="font-size:12pt;font-weight:bold" align="center" | |||
| width="65" height="30" | Perfect match | |||
| width="65" | 1 ins/del | |||
| width="65" | 1 sub | |||
| width="65" | 2 ins/del | |||
| width="65" | 2 sub | |||
| width="65" | 1 ins/del & 1 sub | |||
| align="center" width="65" | 3+ | |||
|- style="font-size:12pt" align="center" valign="bottom" | |||
| align="center" height="15" | 4324 | |||
| align="center" | 2 | |||
| align="center" | 5451 | |||
| align="center" | 0 | |||
| align="center" | 4204 | |||
| align="center" | 22 | |||
| align="center" | 1988 | |||
|} | |||
So it counted WAY fewer probes than were present, and the error rate of substitutions seems to be over estimated (2.48%), where as the error rate of insertions/deletions seem to be underrepresented (0.003% insertions, 1e-4% deletions) | |||
==1% Substitutions== | |||
Generated Data with 1% subsitutions using ''MockHL155_Master, Switch 2''. Copied over to Genome Miner and ran through the bash script from above. | |||
'''Scp''' | |||
scp v4s1_mockseq_error_1s.fq djacobse@132.239.135.41:/media/LTS_15T/DEJ_LTS/SeqStore/130628_HL155/mockseq/error1s/ | |||
===Results=== | |||
2959000 reads; of these: | |||
2959000 (100.00%) were unpaired; of these: | |||
2936721 (99.25%) aligned 0 times | |||
22247 (0.75%) aligned exactly 1 time | |||
32 (0.00%) aligned >1 times | |||
0.75% overall alignment rate | |||
So almost none of the sequences aligned, even with only 1% substitutions. Going to verify that the substitutions worked. |
Latest revision as of 19:28, 4 February 2014
Mock HL155 (Started 9/9/2013)[edit]
Bowtie on Perfect Mock Sequence Data with Perfect Base Qualities[edit]
Based on the results from yesterday I'm seeing what kind of alignment bowtie2 does on perfect sequences with perfect base qualities. Yesterday's sequences were direct copies of the .fa file, but they used base qualities randomly picked from the original sequencing data. These base qualities are all I (an extremely high quality score for phred33).
Generated using MATLAB (MockHL155_Master, Switch 3)
scp to Genome Miner
scp V4S1_mockseq_perfseqbqual.fq djacobse@132.239.135.41:/media/LTS_15T/DEJ_LTS/SeqStore/130628_HL155/mockseq/
bowtie2
bowtie2 -x bbuild_v4s1 -U v4s1_mockseq_perfseqbqual.fq --phred33 -S mockseq_v4s1_perfseqbqual.sam
Alignment Results
2959000 reads; of these: 2959000 (100.00%) were unpaired; of these: 0 (0.00%) aligned 0 times 2953000 (99.80%) aligned exactly 1 time 6000 (0.20%) aligned >1 times 100.00% overall alignment rate
So it looks like even with indisputable base quality and sequences, about 0.2% of sequences are so close that bowtie2 thinks they can be several sequences.
Convert Sam to Bam and Sort
samtools view -bS mockseq_v4s1_perfseqbqual.sam | samtools sort - mockseq_v4s1_perfseqbqual_sort
Filter and Convert Back to Sam
samtools view -h -F 4 -q 8 mockseq_v4s1_perfseqbqual_sort.bam > mockseq_v4s1_perfseqbqual_sortfilt.sam
Replace Matched Pairs with Equals
samtools calmd -eS mockseq_v4s1_perfseqbqual_sortfilt.sam V4S1.fa > mockseqv4s1perfseqbqual_sortfilteq.sam
Perl Script Results script: File:Imp count mismatch.txt
Error Rate: 0 Error Rate of Insertions: 0 Error Rate of Deletions: 0 Error Rate of Substitutions: 0 Count: 2955000 perfect, 0 errors
Bowtie on Perfect Mock Sequence Data (Started 09/11/2013)[edit]
Using Matt's Perlscript: File:Imp count mismatch.txt
Error Rate: 0 Error Rate of Insertions: 0 Error Rate of Deletions: 0 Error Rate of Substitutions: 0 Counts: 2955000 perfect, 0 with error
So this is expected from the sequencing data with no errors. Interestingly, 4000 probes were not counted at all, meaning (most likely) that 4000 probes were of quality too low (they were filtered out in the calmd stage). 4000 is a suspicious number, being even, which could mean that 4 different probes were not counted. It could also just be coincidence. However, the error is present in reads of good quality and reads of normal quality.
Bash Script[edit]
To save myself a lot of time as I run through all of this stuff, I've created a bash script to run through the pipeline (above), going through bowtie2 all the way to the sam file with the equals signs. It's here: Media:Hl155bash.sh
Analysis on Mimic Data[edit]
Mimic data is named such because it has the same error rates calculated from the actual sequencing data (1.01% substitutions, 0.26% insertions, and 0.05% deletions). It was generated using MockHL155_Master, Switch 2 with parameters 1.01, 0.26, and 0.05.
The analysis was done using the bash script from above.
Perl Script Results[edit]
Used Matt's perl script (above)
Error Rate: 0.024866487399162 Error Rate of Insertions: 3.25182915389907e-05 Error Rate of Deletions: 1.25070352073041e-06 Error Rate of Substitutions: 0.0248327184041023
Errors:
Perfect match | 1 ins/del | 1 sub | 2 ins/del | 2 sub | 1 ins/del & 1 sub | 3+ |
4324 | 2 | 5451 | 0 | 4204 | 22 | 1988 |
So it counted WAY fewer probes than were present, and the error rate of substitutions seems to be over estimated (2.48%), where as the error rate of insertions/deletions seem to be underrepresented (0.003% insertions, 1e-4% deletions)
1% Substitutions[edit]
Generated Data with 1% subsitutions using MockHL155_Master, Switch 2. Copied over to Genome Miner and ran through the bash script from above.
Scp
scp v4s1_mockseq_error_1s.fq djacobse@132.239.135.41:/media/LTS_15T/DEJ_LTS/SeqStore/130628_HL155/mockseq/error1s/
Results[edit]
2959000 reads; of these: 2959000 (100.00%) were unpaired; of these: 2936721 (99.25%) aligned 0 times 22247 (0.75%) aligned exactly 1 time 32 (0.00%) aligned >1 times 0.75% overall alignment rate
So almost none of the sequences aligned, even with only 1% substitutions. Going to verify that the substitutions worked.