Editing
Kun:LabNotes/Haplotyping/2011-5-12
(section)
Jump to navigation
Jump to search
Warning:
You are not logged in. Your IP address will be publicly visible if you make any edits. If you
log in
or
create an account
, your edits will be attributed to your username, along with other benefits.
Anti-spam check. Do
not
fill this in!
===Data analysis=== ====Part I==== *Read mapping summary. {| {{table}} | align="center" style="background:#f0f0f0;"|'''Index''' | align="center" style="background:#f0f0f0;"|'''Sample Description''' | align="center" style="background:#f0f0f0;"|'''Total reads''' | align="center" style="background:#f0f0f0;"|'''# Uniquely mapped''' | align="center" style="background:#f0f0f0;"|'''# Failed to map''' | align="center" style="background:#f0f0f0;"|'''# Non-uniquely mapped''' | align="center" style="background:#f0f0f0;"|'''% Uniquely mapped''' | align="center" style="background:#f0f0f0;"|'''% Failed to map''' | align="center" style="background:#f0f0f0;"|'''% Non-uniquely mapped''' | align="center" style="background:#f0f0f0;"|'''# Bis_mapped''' | align="center" style="background:#f0f0f0;"|'''% Bis_mapped''' |- | Indx7||bisulfite-haplotyping-5-cells-1||4,535,639||1,148,608||2,983,032||403,999||25.3%||65.8%||8.9%||2,463,680||54.3% |- | Indx8||bisulfite-haplotyping-5-cells-2||5,908,247||1,350,371||4,062,938||494,938||22.9%||68.8%||8.4%||3,813,863||64.6% |- | Indx9||bisulfite-haplotyping-5-cells-3||4,281,883||689,409||3,333,689||258,785||16.1%||77.9%||6.0%||2,895,651||67.6% |- | Indx10||bisulfite-haplotyping-5-cells-4||3,922,281||953,993||2,630,259||338,029||24.3%||67.1%||8.6%||2,193,243||55.9% |- | Indx11||bisulfite-haplotyping-5-cells-5||3,353,075||451,253||2,737,656||164,166||13.5%||81.6%||4.9%||1,654,575||49.3% |- | Indx12||bisulfite-haplotyping-5-cells-6||4,385,132||678,617||3,449,632||256,883||15.5%||78.7%||5.9%||2,429,874||55.4% |- | |} *On average 27% of the reads can be mapped to the human genome, suggesting that there were contamination of human DNA after bisulfite conversion, most likely at the MDA step or the washing step before MDA. *An average of 58% of the reads can be mapped to the bisulfite converted genome, indicating that roughly 58%-27%=31% of the sequencing reads came from the bisulfite gel trapped chromosomes. Note that genomic reads are mappable to the bisulfite genome because during the mapping the sequencing reads were bisulfite converted in silico. *To further confirm that a significant fraction of sequencing reads are bisulfite converted reads, I did the base composition analysis on the raw reads. I took the top 10,000 reads from the raw fastq file, then calculate the ratio of C/T and A/G for each read. Then I made a scatter plot of C/T versus A/G. I also did the same analysis on a genomic library and a BSPP library for comparison. Clearly C/T or A/G is highly skewed in the bisulfite converted reads, compared with regular genomic reads. The scatter plot indicates that the in-gel methylHap reads is a mixture of regular genomic reads and bisulfite converted reads. [[Image:GenomicReadBaseComposition.png|300px]][[Image:BSPPReadBaseComposition.png|300px]][[Image:InGelMethylHapReadBaseComposition.png|300px]] *This above analysis led to the next question: can we separate bisulfite converted reads from regular genomic reads in the same sequence file based on the base composition? To do this I calculate a metric, which is the absolute log-transformed ratio of [C/T]/[A/G], and ranked 10000 genomic reads and BSPP reads based on the values. Then I created a scatter plot of this metric versus the rank. This plot can provide a good idea on where can I find a cutoff value that can best distinguish the two types of reads. [[Image:BaseBias_genomic_vs_BisReads.png]] Guided by this plot, I found that with a threshold of 2.5, I can keep 98.9% (sensitivity) of the bisulfite reads while 11% (false positives) of the genomic reads will be also included. If I raise the threshold to 3.0, I can still keep 96.4% of bisulfite reads, and the fraction of genomic reads included dropped to 7%. *Alternatively, I can map the read to the human genome first, than perform bisulfite mapping only on the unmappable reads. *To Do: **Extract bisulfite converted reads using the two methods mentioned above, perform bisulfite read mapping. **Compare the methylation level with GM20431 BSPP data. **Calculate bisulfite conversion rate. **Check whether reads from each library tend to concentrate to a subset of genomic regions. If that's the case, it's a good indication that large chromosome fragments were present in different gel slices and therefore we might be able to get haplotype data. ** Perform SNP calling, and compare the genotypes with the PGP1 CGI data. ====Part II==== *Subtraction of genomic reads on all six read files using bowtie: /home/kunzhang/softwares/bowtie-latest/bowtie -k 1 -l 32 -n 2 -p 4--un lane1_Indx7_bisReads.fastq /home/kunzhang/HsGenome/1KG.ref/human_b36_male.fa lane1_Indx7.fastq lane1_Indx7.hg18.bowtie.out *Call methylation on the subtracted bisReads. *I estimated the bisulfite conversion rates by checking the average methylation level on CHH sites. Since CHH methylation is very low, most of the Cs on CHH sites are due to incomplete conversion. [[Media:pileup2ChhMethylLevel.txt|pileup2ChhMethylLevel.pl]] lane1_Indx7_bisReads.fastq.fwd.pileup lane1_Indx7_bisReads.fastq.rev.pileup chr1 {| {{table}} | align="center" style="background:#f0f0f0;"|'''Index''' | align="center" style="background:#f0f0f0;"|'''Total reads''' | align="center" style="background:#f0f0f0;"|'''# reads after genomic filtering''' | align="center" style="background:#f0f0f0;"|'''# reads mapped to bisCvt hg18''' | align="center" style="background:#f0f0f0;"|'''% mapped''' | align="center" style="background:#f0f0f0;"|'''% bisCvt rate''' |- | Indx7||4,535,639||2,990,075||1,119,588||37.44%||78.79% |- | Indx8||5,908,247||4,075,544||2,227,644||54.66%||85.63% |- | Indx9||4,281,883||3,344,978||2,085,038||62.33%||84.94% |- | Indx10||3,922,281||2,635,968||1,066,949||40.48%||76.74% |- | Indx11||3,353,075||2,743,982||1,131,382||41.23%||77.95% |- | Indx12||4,385,132||3,456,344||1,636,709||47.35%||90.23% |} *Comparison of methylation level with BSPP data from four PGP lymphocyte lines. {| {{table}} | align="center" style="background:#f0f0f0;"|'''Correlation''' | align="center" style="background:#f0f0f0;"|'''Indx7''' | align="center" style="background:#f0f0f0;"|'''Indx8''' | align="center" style="background:#f0f0f0;"|'''Indx9''' | align="center" style="background:#f0f0f0;"|'''Indx10''' | align="center" style="background:#f0f0f0;"|'''Indx11''' | align="center" style="background:#f0f0f0;"|'''Indx12''' | align="center" style="background:#f0f0f0;"|'''PGP1''' | align="center" style="background:#f0f0f0;"|'''PGP5''' | align="center" style="background:#f0f0f0;"|'''PGP7''' | align="center" style="background:#f0f0f0;"|'''PGP8''' |- | Indx7||1.000||-||-||-||-||-||-||-||-||- |- | Indx8||0.374||1.000||-||-||-||-||-||-||-||- |- | Indx9||0.637||0.509||1.000||-||-||-||-||-||-||- |- | Indx10||0.035||0.185||0.414||1.000||||-||-||-||-||- |- | Indx11||0.349||0.689||1.000||-0.262||1.000||-||-||-||-||- |- | Indx12||-0.119||-0.465||-0.052||-0.225||-0.211||1.000||-||-||-||- |- | PGP1||0.694||0.511||0.591||0.473||0.526||0.615||1.000||-||-||- |- | PGP5||0.421||0.437||0.708||0.381||0.725||0.562||0.881||1.000||-||- |- | PGP7||0.590||0.368||0.645||0.470||0.715||0.445||0.893||0.904||1.000||- |- | PGP8||0.773||0.508||0.641||0.461||0.590||0.496||0.884||0.932||0.909||1.000 |- | |} *SNP calling. **Raw SNPs were called as part of the methylation mapping process, then filtered. I got between 293-703 SNPs per library. Note that only SNPs with the minimal read depth of 8 could passed the filter. The cut-off for read depth could be reduced in the future. **I wrote a script to compare SNPs called from bisulfite reads to Complete Genomics' data in var files. Roughly 2/3 of the SNPs could be mapped to CGI's data. The majority of such SNPs showed one allele in the bisulfite data, and two alleles in CGI's data, which is expected. A small percentage of SNPs are homozygous in both, but different. When I compared the CGI's genotypes generated from three different cell types of PGP1, it became obvious that these SNPs were called heterozygous in one data set, but homozygous in another, probably due to the lack of read depth or random allelic drops in CGI's data.
Summary:
Please note that all contributions to ZhangLabWiki may be edited, altered, or removed by other contributors. If you do not want your writing to be edited mercilessly, then do not submit it here.
You are also promising us that you wrote this yourself, or copied it from a public domain or similar free resource (see
ZhangLabWiki:Copyrights
for details).
Do not submit copyrighted work without permission!
Cancel
Editing help
(opens in new window)
Navigation menu
Personal tools
Not logged in
Talk
Contributions
Create account
Log in
Namespaces
Page
Discussion
English
Views
Read
Edit
View history
More
Search
Navigation
Main Page
Current events
Recent changes
Random page
Investigators
Matt Cai
Song Chen
Eric Chu
Dinh Diep
Elizabeth Duong
Shicheng Guo
Alan Fung
Daniel Jacobsen
Blue Lake
Huy Lam
Alice Li
Andrew Richards
Brandon Sos
Chris Wei
Yan Wu
Kun Zhang
Tools
What links here
Related changes
Special pages
Page information