Dinh/Dinh 2013/NOTES/2013-10-7: Difference between revisions
Jump to navigation
Jump to search
>Dinh m (→Conclusions) |
>Dinh |
||
(21 intermediate revisions by the same user not shown) | |||
Line 16: | Line 16: | ||
2. Technical replicates for Blueprint samples #1 and #2 from the Aug 28 MiSeq run. (compare different captures) - BP1rep, BP2rep | 2. Technical replicates for Blueprint samples #1 and #2 from the Aug 28 MiSeq run. (compare different captures) - BP1rep, BP2rep | ||
3. Technical replicates for Blueprint samples #1 and #2 from the Sep 13 HiSeq Rapid run. (compare different runs) - BP1rep-hiseq, BP2rep-hiseq | 3. Technical replicates for Blueprint samples #1 and #2 from the Sep 13 HiSeq Rapid run. (compare different runs) - BP1rep-hiseq, BP2rep-hiseq | ||
4. NA12878 captured with the original probes set and sequenced in HL157 (+ control) - HMrep1 | |||
5. NA12878 captured with the normalized probes set and sequenced in HL157 (+ control) - HMrep2 | 5. NA12878 captured with the normalized probes set and sequenced in HL157 (+ control) - HMrep2 | ||
6. NA12878 captured with the normalized probes set and sequenced in HiSeq Rapid run - HMrep3 | 6. NA12878 captured with the normalized probes set and sequenced in HiSeq Rapid run - HMrep3 | ||
7. DMR330K capture from the Mar22 MiSeq run as positive control. (+ control) - H1-dmr330k | 7. DMR330K capture from the Mar22 MiSeq run as positive control. (+ control) - H1-dmr330k | ||
* All reads were sampled from the first | 8. NA12878 captured with the original probes set and sequenced in June 20 MiSeq run (+ control) - HM-Miseq1, HM-Miseq2 | ||
* All reads were sampled from only read 1. | |||
* Sampling was done randomly by using "sort -R" and then taking the first 1000 reads after random sorting. | |||
* All reads were trimmed 27 bp from the 5' end and quality trimmed using mean quality score of 30. | * All reads were trimmed 27 bp from the 5' end and quality trimmed using mean quality score of 30. | ||
Line 27: | Line 29: | ||
* Library free adapters/linkers/primers | * Library free adapters/linkers/primers | ||
>AmpF6.3Sol | >AmpF6.3Sol | ||
AATGATACGGCGACCACCGACACTCT<b>CAGATGTTATCGAGGTCCGAC</b> | |||
>AmpR6.3Ind1 | >AmpR6.3Ind1 | ||
CAAGCAGAAGACGGCATACGAGATCGTGATGC<b>TAGGAACGATGAGCCTCCAA</b>C | |||
>SolSeqV6.3.2 | >SolSeqV6.3.2 | ||
ACCGACACTCT<b>CAGATGTTATCGAGGTCCGAC</b> | |||
>SolSeqV6.3.2r | >SolSeqV6.3.2r | ||
GC<b>TAGGAACGATGAGCCTCCAAC</b> | |||
>AmpR6.3IndSeq | >AmpR6.3IndSeq | ||
<b>GTTGGAGGCTCATCGTTCCTA</b>GC | |||
>AmpF6.4Sol | >AmpF6.4Sol | ||
AATGATACGGCGACCACCGAGATCTACACCACTCT<b>CAGATGTTATCGAGGTCCGAC</b> | |||
>SolSeqV6.3.3 | >SolSeqV6.3.3 | ||
TACACCACTCT<b>CAGATGTTATCGAGGTCCGAC</b> | |||
* BLAST args: | * BLAST args: | ||
$blastn -p blastn -d $db_file -K 100 -e 0.0001 -g F -F "m D" -a 4 -W 7 -m 8 -i $f -o $f.primers.blast | $blastn -p blastn -d $db_file -K 100 -e 0.0001 -g F -F "m D" -a 4 -W 7 -m 8 -i $f -o $f.primers.blast | ||
Line 50: | Line 52: | ||
222 BP1rep-hiseq | 222 BP1rep-hiseq | ||
558 BP2rep-hiseq | 558 BP2rep-hiseq | ||
'''12''' HMrep1 (but this library has ~18.70% mapping rate??) | |||
'''5''' HMrep2 (but this library has ~15.07% mapping rate??) | '''5''' HMrep2 (but this library has ~15.07% mapping rate??) | ||
505 HMrep3-hiseq | 505 HMrep3-hiseq | ||
'''23''' H1-dmr330K | '''23''' H1-dmr330K | ||
* Number of hits to each primer sequences. NOTE: a read can hit multiple primer sequences. | 279 HM-Miseq1 | ||
332 HM-Miseq2 | |||
* Number of hits to each primer sequences. NOTE: a read can hit multiple primer sequences. In parentheses are regions within primers which was hit 99% of the time. | |||
{| {{table}} border=1 | {| {{table}} border=1 | ||
| align="center" style="background:#f0f0f0;"| | | align="center" style="background:#f0f0f0;"|'''Sequencing Length''' | ||
| align="center" style="background:#f0f0f0;"|'''Sample''' | |||
| align="center" style="background:#f0f0f0;"|'''AmpF6.3Sol''' | | align="center" style="background:#f0f0f0;"|'''AmpF6.3Sol''' | ||
| align="center" style="background:#f0f0f0;"|'''AmpR6.3Ind1''' | | align="center" style="background:#f0f0f0;"|'''AmpR6.3Ind1''' | ||
Line 65: | Line 70: | ||
| align="center" style="background:#f0f0f0;"|'''SolSeqV6.3.3''' | | align="center" style="background:#f0f0f0;"|'''SolSeqV6.3.3''' | ||
|- | |- | ||
| BP1||510||555||590||555||555||510||510 | | 250 || BP1||510 (27->47)||555 (53->33)||590 (12->32)||555 (23->3)||555 (1->21)||510 (36->56)||510 (12->32) | ||
|- | |- | ||
| BP2||581||639||667||639||639||581||581 | | 250 || BP2||581||639||667||639||639||581||581 | ||
|- | |- | ||
| BP1rep||363||444||458||444||444||363||363 | | 250 || BP1rep||363||444||458||444||444||363||363 | ||
|- | |- | ||
| BP2rep||653||786||811||785||785||653||653 | | 250 || BP2rep||653||786||811||785||785||653||653 | ||
|- | |- | ||
| | | 250 || HM-Miseq1 ||245||270||278||270||270||245||245 | ||
|- | |- | ||
| | | 250 || HM-Miseq2 ||289||321||331||320||320||289||289 | ||
|- | |- | ||
| | | 100 || BP1rep-hiseq||8||220||222||220||220||8||8 | ||
|- | |- | ||
| | | 100 || BP2rep-hiseq||12||556||558||556||556||12||12 | ||
|- | |- | ||
| | | 100 || HMrep3-hiseq||10||502||505||502||502||10||10 | ||
|- | |- | ||
| | | 150 || H1-dmr330K||12||23||23||23||23||12||12 | ||
|} | |- | ||
| 60 || HMrep1||0||5||5||5||5||0||0 | |||
|- | |||
| 60 || HMrep2||1||12||12||12||12||1||1 | |||
|} | |||
* The reads map to ONLY the adapter sequences. | |||
* Frequency of positions where the adapter mapped to reads: | |||
[[File:blueprint.blast.primers_positions.png]] | |||
===Conclusions=== | ===Conclusions=== | ||
* Chimeric sequences were present in all except for the H1-dmr330k | * Chimeric sequences were present in all except for the H1-dmr330k which was 124bp long after all the necessary trimmings. The samples sequenced in HL157 was only 35bp after trimming, which might affected the detection of the primer sequences. Since in the Hiseq run with 65bp, we can already start to see the the primer sequences, which was not present in the H1-dmr330k sample, there is significant mapping of reads to the adapter/linker/primer sequences with the Blueprint probes compared to the DMR330k probes. | ||
== BLAST sequences to NT database == | == BLAST sequences to NT database == | ||
* If capture reagents were contaminated, we would expect that reads from capture of the converted sample would not map to the unconverted NT database. Reads that map to the unconverted NT database would means that we have contaminated capture reagents. | * If capture reagents were contaminated, we would expect that reads from capture of the converted sample would not map to the unconverted NT database. Reads that map to the unconverted NT database would means that we have contaminated capture reagents. | ||
* | * BLAST args: | ||
$blastn -p blastn -d $dbfile -K 100 -e 100 -g F -F "m D" -a 4 -W 7 -i $f.fa -o $f.nt.blast | |||
* Results: No evident of significant contamination from DMR330K capture with H1 sample. Same for BP capture with HapMap and Blueprint samples. | |||
{| {{table}} border=1 | |||
| align="center" style="background:#f0f0f0;"|'''Sample''' | |||
| align="center" style="background:#f0f0f0;"|'''Graph''' | |||
| align="center" style="background:#f0f0f0;"|'''No hits''' | |||
| align="center" style="background:#f0f0f0;"|'''Not assigned''' | |||
| align="center" style="background:#f0f0f0;"|'''Hominidae''' | |||
|- | |||
| BP1||[[File:MEGAN.BP.mapNT.BP1.png | 200px ]]||17||482||347 | |||
|- | |||
| BP2||[[File:MEGAN.BP.mapNT.BP2.png | 200px ]]||17||533||307 | |||
|- | |||
| BP1rep||[[File:MEGAN.BP.mapNT.BP1rep.png | 200px ]]||25||468||309 | |||
|- | |||
| BP2rep||[[File:MEGAN.BP.mapNT.BP2rep.png | 200px ]]||17||573||197 | |||
|- | |||
| H1-dmr330K||[[File:MEGAN.BP.mapNT.H1-dmr330k.png | 200px ]]||54||422||356 | |||
|- | |||
| BP1rep-hiseq||[[File:MEGAN.BP.mapNT.BP1rep-hiseq.png | 200px ]]||138||576||169 | |||
|- | |||
| BP2rep-hiseq||[[File:MEGAN.BP.mapNT.BP2rep-hiseq.png | 200px ]]||53||758||91 | |||
|- | |||
| Hmrep3-hiseq||[[File:MEGAN.BP.mapNT.HMrep3.png | 200px ]]||55||724||110 | |||
|- | |||
| Hmrep1||[[File:MEGAN.BP.mapNT.HMrep1.png | 200px ]]||42||868||24 | |||
|- | |||
| Hmrep2||[[File:MEGAN.BP.mapNT.HMrep2.png | 200px ]]||45||852||44 | |||
|- | |||
| Hm-Miseq||[[File:MEGAN.BP.mapNT.HM-miseq.png | 200px ]]||29||385||438 | |||
|} | |||
== BLAST sequences to RepBase == | == BLAST sequences to RepBase == | ||
* If there were off-target capture to repeats, we would get hits to the converted RepBase sequences, but they would be mostly human. | * If there were off-target capture to repeats, we would get hits to the converted RepBase sequences, but they would be mostly human. | ||
Line 124: | Line 168: | ||
| H1-dmr330K*||980||1||20||549||33||449 | | H1-dmr330K*||980||1||20||549||33||449 | ||
|- | |- | ||
| | | HMrep2||990||4||10||702||37||287 | ||
|- | |- | ||
| | | HMrep1||991||1||7||750||35||249 | ||
|- | |- | ||
| BP1rep-hiseq||972||16||25||640||78||354 | | BP1rep-hiseq||972||16||25||640||78||354 | ||
Line 132: | Line 176: | ||
| BP2rep-hiseq||976||14||17||660||51||334 | | BP2rep-hiseq||976||14||17||660||51||334 | ||
|- | |- | ||
| | | HMrep3||976||15||21||652||71||342 | ||
|- | |||
| HM-Miseq1||967||12||24||458||57||540 | |||
|- | |||
| HM-Miseq2||976||5||22||471||67||524 | |||
|} | |} | ||
Latest revision as of 17:55, 10 October 2013
Blueprint Capture, low mapping issue[edit]
- Low mapping can be caused by any of the following problems
- Contamination - contamination of library after capture would not be amplified by the last PCR step or sequenced due to lack of adapter sequences. Contamination of reagents during capture may lead to low mapping rates. Contamination may explain why we have such variable mapping rates, since each DNA sample had varying quality, it's possible that the lower quality samples would have more capture from contaminating sequences and thus lower mapping rate.
- Do the amount of input converted gDNA correlate with mapping rate?
- Do the captured sequences map to other genomes in NT database?
- Chimeric capture sequences - When imperfect hybridization of capture arms to genomic sequences which occurs frequently in the genome, leading to off-target & chimeric sequences. Or when the capture arms capture each other.
- How many of the sequences map to the adapters, linkers, and primers?
- Do the sequences map to human repeats? - No, when mapping with Bowtie, most of the reads either map uniquely or not at all.
- Do the sequences map to non-human repeats? - Check by performing Blast to RepBase18.08
- Error from sequencing - We did not use a PhiX control for the MiSeq runs, but had a control for the HiSeq RapidRun (control was ~ 89% of flowcell, BP libraries was ~ 11%). The mapping rate of the HiSeq RapidRun was still low for many samples. Also, the majority of mapped reads have low error rates at the 3' end.
- Error from oligo synthesis - error from oligo synthesis would still allow us to recover the alignments from the off-target regions. Even so, we should get lower the specificity of capture, but this is not the case.
Reads sampling and trimming[edit]
- Randomly sample 1000 reads from the following libraries for analysis:
1. Blueprint samples #1 and #2 from the Aug 28 MiSeq run. (compare different BP samples) - BP1, BP2 2. Technical replicates for Blueprint samples #1 and #2 from the Aug 28 MiSeq run. (compare different captures) - BP1rep, BP2rep 3. Technical replicates for Blueprint samples #1 and #2 from the Sep 13 HiSeq Rapid run. (compare different runs) - BP1rep-hiseq, BP2rep-hiseq 4. NA12878 captured with the original probes set and sequenced in HL157 (+ control) - HMrep1 5. NA12878 captured with the normalized probes set and sequenced in HL157 (+ control) - HMrep2 6. NA12878 captured with the normalized probes set and sequenced in HiSeq Rapid run - HMrep3 7. DMR330K capture from the Mar22 MiSeq run as positive control. (+ control) - H1-dmr330k 8. NA12878 captured with the original probes set and sequenced in June 20 MiSeq run (+ control) - HM-Miseq1, HM-Miseq2
- All reads were sampled from only read 1.
- Sampling was done randomly by using "sort -R" and then taking the first 1000 reads after random sorting.
- All reads were trimmed 27 bp from the 5' end and quality trimmed using mean quality score of 30.
BLAST sequences to adapters, linkers, and primers[edit]
- Since these were not bisulfite converted, I use the original reads to blast.
- Library free adapters/linkers/primers
>AmpF6.3Sol AATGATACGGCGACCACCGACACTCTCAGATGTTATCGAGGTCCGAC >AmpR6.3Ind1 CAAGCAGAAGACGGCATACGAGATCGTGATGCTAGGAACGATGAGCCTCCAAC >SolSeqV6.3.2 ACCGACACTCTCAGATGTTATCGAGGTCCGAC >SolSeqV6.3.2r GCTAGGAACGATGAGCCTCCAAC >AmpR6.3IndSeq GTTGGAGGCTCATCGTTCCTAGC >AmpF6.4Sol AATGATACGGCGACCACCGAGATCTACACCACTCTCAGATGTTATCGAGGTCCGAC >SolSeqV6.3.3 TACACCACTCTCAGATGTTATCGAGGTCCGAC
- BLAST args:
$blastn -p blastn -d $db_file -K 100 -e 0.0001 -g F -F "m D" -a 4 -W 7 -m 8 -i $f -o $f.primers.blast
- Number of reads mapped to library_free_primers:
(BP - Blueprint sample, HM - HapMap sample NA12878, H1 - Human ESC H1) 590 BP1 667 BP2 458 BP1rep 812 BP2rep 222 BP1rep-hiseq 558 BP2rep-hiseq 12 HMrep1 (but this library has ~18.70% mapping rate??) 5 HMrep2 (but this library has ~15.07% mapping rate??) 505 HMrep3-hiseq 23 H1-dmr330K 279 HM-Miseq1 332 HM-Miseq2
- Number of hits to each primer sequences. NOTE: a read can hit multiple primer sequences. In parentheses are regions within primers which was hit 99% of the time.
Sequencing Length | Sample | AmpF6.3Sol | AmpR6.3Ind1 | SolSeqV6.3.2 | SolSeqV6.3.2r | AmpR6.3IndSeq | AmpF6.4Sol | SolSeqV6.3.3 |
250 | BP1 | 510 (27->47) | 555 (53->33) | 590 (12->32) | 555 (23->3) | 555 (1->21) | 510 (36->56) | 510 (12->32) |
250 | BP2 | 581 | 639 | 667 | 639 | 639 | 581 | 581 |
250 | BP1rep | 363 | 444 | 458 | 444 | 444 | 363 | 363 |
250 | BP2rep | 653 | 786 | 811 | 785 | 785 | 653 | 653 |
250 | HM-Miseq1 | 245 | 270 | 278 | 270 | 270 | 245 | 245 |
250 | HM-Miseq2 | 289 | 321 | 331 | 320 | 320 | 289 | 289 |
100 | BP1rep-hiseq | 8 | 220 | 222 | 220 | 220 | 8 | 8 |
100 | BP2rep-hiseq | 12 | 556 | 558 | 556 | 556 | 12 | 12 |
100 | HMrep3-hiseq | 10 | 502 | 505 | 502 | 502 | 10 | 10 |
150 | H1-dmr330K | 12 | 23 | 23 | 23 | 23 | 12 | 12 |
60 | HMrep1 | 0 | 5 | 5 | 5 | 5 | 0 | 0 |
60 | HMrep2 | 1 | 12 | 12 | 12 | 12 | 1 | 1 |
- The reads map to ONLY the adapter sequences.
- Frequency of positions where the adapter mapped to reads:
File:Blueprint.blast.primers positions.png
Conclusions[edit]
- Chimeric sequences were present in all except for the H1-dmr330k which was 124bp long after all the necessary trimmings. The samples sequenced in HL157 was only 35bp after trimming, which might affected the detection of the primer sequences. Since in the Hiseq run with 65bp, we can already start to see the the primer sequences, which was not present in the H1-dmr330k sample, there is significant mapping of reads to the adapter/linker/primer sequences with the Blueprint probes compared to the DMR330k probes.
BLAST sequences to NT database[edit]
- If capture reagents were contaminated, we would expect that reads from capture of the converted sample would not map to the unconverted NT database. Reads that map to the unconverted NT database would means that we have contaminated capture reagents.
- BLAST args:
$blastn -p blastn -d $dbfile -K 100 -e 100 -g F -F "m D" -a 4 -W 7 -i $f.fa -o $f.nt.blast
- Results: No evident of significant contamination from DMR330K capture with H1 sample. Same for BP capture with HapMap and Blueprint samples.
Sample | Graph | No hits | Not assigned | Hominidae |
BP1 | File:MEGAN.BP.mapNT.BP1.png | 17 | 482 | 347 |
BP2 | File:MEGAN.BP.mapNT.BP2.png | 17 | 533 | 307 |
BP1rep | File:MEGAN.BP.mapNT.BP1rep.png | 25 | 468 | 309 |
BP2rep | File:MEGAN.BP.mapNT.BP2rep.png | 17 | 573 | 197 |
H1-dmr330K | File:MEGAN.BP.mapNT.H1-dmr330k.png | 54 | 422 | 356 |
BP1rep-hiseq | File:MEGAN.BP.mapNT.BP1rep-hiseq.png | 138 | 576 | 169 |
BP2rep-hiseq | File:MEGAN.BP.mapNT.BP2rep-hiseq.png | 53 | 758 | 91 |
Hmrep3-hiseq | File:MEGAN.BP.mapNT.HMrep3.png | 55 | 724 | 110 |
Hmrep1 | File:MEGAN.BP.mapNT.HMrep1.png | 42 | 868 | 24 |
Hmrep2 | File:MEGAN.BP.mapNT.HMrep2.png | 45 | 852 | 44 |
Hm-Miseq | File:MEGAN.BP.mapNT.HM-miseq.png | 29 | 385 | 438 |
BLAST sequences to RepBase[edit]
- If there were off-target capture to repeats, we would get hits to the converted RepBase sequences, but they would be mostly human.
- Reads were encoded by converting G to A (all read 1 would be on the complement strand)
- BLAST args:
$blastn -p blastn -d $db_dir/all_repeats.bis.fwd -K 100 -e 0.0001 -g F -F "m D" -a 4 -W 7 -i $f.encoded_2 -o $f.bis.fwd.all_repeats.blast
' | Unconverted repeats | Unconverted repeats | Unconverted repeats | Bis-forward strand | Bis-forward strand | Bis-forward strand |
Unmapped reads | Human repeat reads | Non-human repeat reads | Unmapped reads | Human repeat reads | Non-human repeat reads | |
BP1 | 977 | 7 | 19 | 396 | 94 | 602 |
BP2 | 970 | 11 | 29 | 438 | 89 | 558 |
BP1rep | 930 | 35 | 59 | 455 | 131 | 538 |
BP2rep | 967 | 16 | 24 | 441 | 90 | 553 |
H1-dmr330K* | 980 | 1 | 20 | 549 | 33 | 449 |
HMrep2 | 990 | 4 | 10 | 702 | 37 | 287 |
HMrep1 | 991 | 1 | 7 | 750 | 35 | 249 |
BP1rep-hiseq | 972 | 16 | 25 | 640 | 78 | 354 |
BP2rep-hiseq | 976 | 14 | 17 | 660 | 51 | 334 |
HMrep3 | 976 | 15 | 21 | 652 | 71 | 342 |
HM-Miseq1 | 967 | 12 | 24 | 458 | 57 | 540 |
HM-Miseq2 | 976 | 5 | 22 | 471 | 67 | 524 |
- NOTE: A read may hit both human and non-human repeats
Conclusions[edit]
- BP probes tend to have more human repetitive sequences, but the amount of reads mapping to non-human converted/unconverted repetitive sequences appear to be comparable to the DMR330K probes.
- We do not have off-target capture to repetitive sequences.