Dinh/Dinh 2012/NOTES/2012-2-27
Jump to navigation
Jump to search
Mapping bias to reference[edit]
- I extracted all heterozygous SNP calls from: /home/kunzhang/AfricanGenomes/Tishkoff15/dbSNPAnnotated-GS00319-DNA_C02_1100_37-ASM.tsv
awk '{if($10 ~ /het/)print $0}' /home/kunzhang/AfricanGenomes/Tishkoff15/dbSNPAnnotated-GS00319-DNA_C02_1100_37-ASM.tsv > tmp &
- I used this script to prepare chromosome 10 for bisulfite reads mapping:
File:PrepareAltGenome.txt ./prepareAltGenome.pl /media/1TB_store1/Rui_Scratch/iGenome/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa CAPB046_D02.hg19.hetSNP convert=yes context=cg ./prepareAltGenome.pl /media/1TB_store1/Rui_Scratch/iGenome/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa convert=yes context=cg /usr/local/soap/2bwt-builder var.genome.fa.bis.fwd /usr/local/soap/2bwt-builder var.genome.fa.bis.rev /usr/local/soap/2bwt-builder genome.fa.bis.fwd /usr/local/soap/2bwt-builder genome.fa.bis.rev /usr/local/samtools/samtools faidx var.genome.fa.bis.fwd /usr/local/samtools/samtools faidx var.genome.fa.bis.rev /usr/local/samtools/samtools faidx genome.fa.bis.fwd /usr/local/samtools/samtools faidx genome.fa.bis.rev
- Then I mapped using bisReadMapper.pl (using the reads from PennAfrican_Batch1_Indx2_1.txt, note this is only half of the reads)
nohup bisReadMapper_v1/bisReadMapper.pl CAPB46_D02.params > CAPB46_ref.log & nohup ../bisReadMapper_v1/bisReadMapper.pl CAPB46_D02.paramsVAR > CAPB46_var.log &
- BisReadMapper.pl kept all the bam files, so I used samtools to convert to sam and counted the lines:
420926 ref.fwd.sam 416280 ref.rev.sam ================== 421103 var.fwd.sam (0.4% more reads mapped) 454589 var.rev.sam (9.2% more reads mapped)
- Do reads which we missed when mapping to the reverse strand going to affect the final CpG measurements?
- Correlation HuRef CpGs versus CAPB046 CpGs = 0.9997136
bedtools intersect -wao -a ref.methyl.cg-pos.BED.txt -b var.methyl.cg-pos.BED.txt | awk '{if($10 ~ /chr/) print $4"\t"$13;}' > refvsvar.values Calculate correlation using R
- I also compared to Noi's B01_CAPB046.BED.txt file (this used all of the reads.) Just to validate that the new genome preparation script was not different from the pipeline that Noi is using.
- My HuRef CpGs versus Noi's HuRef CpGs = 0.9911043
- Noi's HuRef CpGs versus CAPB046 CpGs = 0.9910447
Conclusion[edit]
- Reference genome mapping bias is only causing ~9% fewer reads getting mapped to the reverse complementary strand.
- Since we convert all Cs to Ts in the forward strand, the reduce complexity genome is less affected by C/T SNPs (since any C/T snps will be treated the same)
- However, for the reverse complementary strand, where there is a C/T SNP, there is a G (mismatch) to the actual A allele.
- Correlation of CpG methylation values indicates little technical bias.