Dinh/Dinh 2013/NOTES/2013-10-7: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Dinh
>Dinh
 
(29 intermediate revisions by the same user not shown)
Line 18: Line 18:
   4. NA12878 captured with the original probes set and sequenced in HL157 (+ control) - HMrep1
   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 read from each run.
  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.
== BLAST sequences to adapters, linkers, and primers==
== BLAST sequences to adapters, linkers, and primers==
* Since these were not bisulfite converted, I use the original reads to blast.
* Since these were not bisulfite converted, I use the original reads to blast.
* Library free adapters/linkers/primers
* Library free adapters/linkers/primers
  >AmpF6.3Sol
  >AmpF6.3Sol
  AATGATACGGCGACCACCGACACTCTCAGATGTTATCGAGGTCCGAC
  AATGATACGGCGACCACCGACACTCT<b>CAGATGTTATCGAGGTCCGAC</b>
  >AmpR6.3Ind1
  >AmpR6.3Ind1
  CAAGCAGAAGACGGCATACGAGATCGTGATGCTAGGAACGATGAGCCTCCAAC
  CAAGCAGAAGACGGCATACGAGATCGTGATGC<b>TAGGAACGATGAGCCTCCAA</b>C
  >SolSeqV6.3.2
  >SolSeqV6.3.2
  ACCGACACTCTCAGATGTTATCGAGGTCCGAC
  ACCGACACTCT<b>CAGATGTTATCGAGGTCCGAC</b>
  >SolSeqV6.3.2r
  >SolSeqV6.3.2r
  GCTAGGAACGATGAGCCTCCAAC
  GC<b>TAGGAACGATGAGCCTCCAAC</b>
  >AmpR6.3IndSeq
  >AmpR6.3IndSeq
  GTTGGAGGCTCATCGTTCCTAGC
  <b>GTTGGAGGCTCATCGTTCCTA</b>GC
  >AmpF6.4Sol
  >AmpF6.4Sol
  AATGATACGGCGACCACCGAGATCTACACCACTCTCAGATGTTATCGAGGTCCGAC
  AATGATACGGCGACCACCGAGATCTACACCACTCT<b>CAGATGTTATCGAGGTCCGAC</b>
  >SolSeqV6.3.3
  >SolSeqV6.3.3
  TACACCACTCTCAGATGTTATCGAGGTCCGAC
  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 53: Line 56:
   505 HMrep3-hiseq
   505 HMrep3-hiseq
   '''23'''  H1-dmr330K
   '''23'''  H1-dmr330K
* Number of hits to each 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 64: 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)
|-
| 250 || BP2||581||639||667||639||639||581||581
|-
| 250 || BP1rep||363||444||458||444||444||363||363
|-
|-
| BP2||581||639||667||639||639||581||581
| 250 || BP2rep||653||786||811||785||785||653||653
|-
|-
| BP1rep||363||444||458||444||444||363||363
| 250 || HM-Miseq1 ||245||270||278||270||270||245||245
|-
|-
| BP2rep||653||786||811||785||785||653||653
| 250 || HM-Miseq2 ||289||321||331||320||320||289||289
|-
|-
| BP1rep-hiseq||8||220||222||220||220||8||8
| 100 || BP1rep-hiseq||8||220||222||220||220||8||8
|-
|-
| BP2rep-hiseq||12||556||558||556||556||12||12
| 100 || BP2rep-hiseq||12||556||558||556||556||12||12
|-
|-
| H1-dmr330K||12||23||23||23||23||12||12
| 100 || HMrep3-hiseq||10||502||505||502||502||10||10
|-
|-
| HMrep1||0||5||5||5||5||0||0
| 150 || H1-dmr330K||12||23||23||23||23||12||12
|-
|-
| HMrep2||1||12||12||12||12||1||1
| 60 || HMrep1||0||5||5||5||5||0||0
|-
|-
| HMrep3-hiseq||10||502||505||502||502||10||10
| 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 expect for the H1-dmr330k run and the HapMap (HM) samples sequenced in HL157.
* 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.
* Waiting for mapping results....
* 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 122: 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
| HMrep2||990||4||10||702||37||287
|-
|-
| Hmrep1*||991||1||7||750||35||249
| HMrep1||991||1||7||750||35||249
|-
|-
| BP1rep-hiseq||972||16||25||640||78||354
| BP1rep-hiseq||972||16||25||640||78||354
Line 130: 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
| 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
* NOTE: A read may hit both human and non-human repeats
===Conclusions===
===Conclusions===
* DMR330K and BP capture in HL157 on HapMap sample have very little of human repeats. All normalized BP capture high % of human repeats.
* 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.

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
  1. 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.
    1. Do the amount of input converted gDNA correlate with mapping rate?
    2. Do the captured sequences map to other genomes in NT database?
  2. 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.
    1. How many of the sequences map to the adapters, linkers, and primers?
    2. Do the sequences map to human repeats? - No, when mapping with Bowtie, most of the reads either map uniquely or not at all.
    3. Do the sequences map to non-human repeats? - Check by performing Blast to RepBase18.08
  3. 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.
  4. 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.