Dinh/Dinh 2012/NOTES/2012-2-27

From ZhangLabWiki
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.