Mutation Rate and Spectra

[hollygene]

Calculate mutation rate and find spectrum of mutations for each strain in dataset.

Begin with H0

[hollygene]

Methods - Following GATK's Best Practices https://github.com/hollygene/TE_MA/blob/TE_MA/H0_script.sh

• Create unmapped bam from paired-end fastqs using picard

• mark illumina adapters

• convert unmapped bams back to fastq

• align to reference using bwa mem

• merge bams using MergeBamAlignment

• sort bam output

• index bam output

• Mark and remove duplicates using picard

• call mutations using HaplotypeCaller in gVCF mode

• combine 8 random samples together and jointly genotype

• use resulting .vcf file to recalibrate quality scores

• re-call variants using haplotype caller with recalibrated samples

• combine gVCFs from each sample

• jointly genotype gVCFs

• filter final VCF

  • remove sites with low depth (double-check the method here)
  • remove sites with high depth
  • remove low mappability sites found from GEM
  • select only those sites that are called in the ancestor (i.e. no ./.)
  • remove sites that are heterozygous in the ancestor
  • keep only sites that are biallelic

• Dave modified the calls in excel using the following codes:

• pull out ancestor calls -MID(M2,SEARCH(":",M2)+1,SEARCH(":",M2,SEARCH(":",M2)+1)-SEARCH(":",M2)-1)

• verify number of ancestor alleles

  • 1+LEN(DD2)-LEN(SUBSTITUTE(DD2,",",""))

• isolate first allele depth

  • IF(DE2=2,LEFT(DD2,SEARCH(",",DD2)-1),"na")

• isolate second allele depth

  • IF(DE2=2,RIGHT(DD2,LEN(DD2)-SEARCH(",",DD2)),"na")

• get depth total

  • IF(DF2="na","",DF2+DG2)

• verify depth between 50 and 170?

  • IF(DH2="","",IF(DH2<50,"low",IF(DH2>170,"high","ok")))

• find Ancestor allele freq

  • IF(DH2="","",DF2/DH2)

• 0.1<ancestor_freq<0.9

  • IF(DJ2="","",IF(AND(DJ2>0.1,DJ2<0.9),"discard","keep"))

• 0.05<ancestor_freq<0.95

  • IF(DJ2="","",IF(AND(DJ2>0.05,DJ2<0.95),"discard","keep"))

• number of "./." in any line or ancestor

  • COUNTIF(BI2:DC2,"./.")

• overall (0.1, 0.9)

  • IF(DM2=0,IF(AND(DI2="ok",DK2="keep"),"keep","discard"),"discard")

• overall (0.05, 0.95

  • IF(DM2=0,IF(AND(DI2="ok",DL2="keep"),"keep","discard"),"discard")

• heterozygous ancestor

  • COUNTIF(BI2,"0/1")

• keep?

  • IF(AND(DM2=0,DP2=0),"keep","discard")

• check (het anc are being removed)(

  • IF(AND(DO2="keep",DQ2="discard"),"problem","ok")

• Ancestor genotype

  • =BI2

• number of hets in MA lines

  • =COUNTIF(BJ2:DC2,"0/1")

• 0/0 in MA lines

  • =COUNTIF(BJ2:DC2,"0/0")

• 1/1 in MA lines

  • =COUNTIF(BJ2:DC2,"1/1")