Dinh/Dinh 2012/NOTES/2012-9-22

From ZhangLabWiki
Revision as of 23:58, 24 September 2012 by >Dinh
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to navigation Jump to search

Use BisSNP for calling SNPs for UCLA SZ datasets[edit]

  • bisReadMapper to map -> bam files -> combine Watson and Crick -> use samtools fillmd -> processed bam file
  • Merge bam files from 4 different batches using samtools, then index and sort by chr positions
  • Run BisSNP on sorted, merged and combined bam files (96 samples in total)
 I ran some of the samples on meangenemachine and genemapster.
 I moved to triton cluster to save time (but need to figure out approximately how much compute time and memory is needed per sample)
 Total time for genotyping is 36 hours per sample. I skipped base recalibration (~ 24 hours) and indel realignment.
  • Make sure that tabix and bgzip2 are in /usr/bin (required for vcf-merge):
 http://samtools.sourceforge.net/tabix.shtml
  • Run perl script, the input file is a sample_info file :
#!/usr/bin/perl -w
use strict;
my $line = <STDIN>;
while(my $line = <STDIN>){
        chomp($line);
        my ($indx, $o_sn, $c_sn, $tfam, $order) = split /\t/, $line;
        # save the header; remove the problematic header line with vcfProcessLog
        system("grep \"^#\" $indx.fillmd.snp.raw.vcf |  awk '{if(\$0 !~ /vcfProcessLog/) print \$0;}' | sed 's/hs/$tfam/g' > $indx.fillmd.snp.raw.vcf.processed");
        # sort the rest, filter by minimum phred score 20
        system("grep -v \"^#\" $indx.fillmd.snp.raw.vcf | ../filterLowQual.pl | sort -k1,1 -k2,2n >> $indx.fillmd.snp.raw.vcf.processed");
        # compress 
        system("bgzip $indx.fillmd.snp.raw.vcf.processed");
        # make tabix index
        system("tabix -p vcf $indx.fillmd.snp.raw.vcf.processed.gz");
        # print names for vcf-merge
        print $indx, ".fillmd.snp.raw.vcf.processed.gz ";
}
  • Run vcftools to merge individual gz files into an all-individual vcf file:
 ~ddiep/softwares/vcftools_0.1.9/bin/vcf-merge <printed list of *.gz files in order> > Merged_VCF
  • Sample names used:
Using column name 'GK000501' for Indx1.fillmd.snp.raw.vcf.processed.gz:GK000501
Using column name 'GK000502' for Indx2.fillmd.snp.raw.vcf.processed.gz:GK000502
Using column name 'GK000503a' for Indx3.fillmd.snp.raw.vcf.processed.gz:GK000503a
Using column name 'GK000504' for Indx4.fillmd.snp.raw.vcf.processed.gz:GK000504
Using column name 'GK002701a' for Indx5.fillmd.snp.raw.vcf.processed.gz:GK002701a
Using column name 'GK002702' for Indx6.fillmd.snp.raw.vcf.processed.gz:GK002702
Using column name 'GK002703' for Indx7.fillmd.snp.raw.vcf.processed.gz:GK002703
Using column name 'GK002704' for Indx8.fillmd.snp.raw.vcf.processed.gz:GK002704
Using column name 'GK005001' for Indx9.fillmd.snp.raw.vcf.processed.gz:GK005001
Using column name 'GK005002' for Indx10.fillmd.snp.raw.vcf.processed.gz:GK005002
Using column name 'GK005003' for Indx11.fillmd.snp.raw.vcf.processed.gz:GK005003
Using column name 'GK005004' for Indx12.fillmd.snp.raw.vcf.processed.gz:GK005004
Using column name 'GK027702' for Indx13.fillmd.snp.raw.vcf.processed.gz:GK027702
Using column name 'GK009101a' for Indx14.fillmd.snp.raw.vcf.processed.gz:GK009101a
Using column name 'GK009102' for Indx15.fillmd.snp.raw.vcf.processed.gz:GK009102
Using column name 'GK009103' for Indx16.fillmd.snp.raw.vcf.processed.gz:GK009103
Using column name 'GK009104' for Indx17.fillmd.snp.raw.vcf.processed.gz:GK009104
Using column name 'GK010001' for Indx18.fillmd.snp.raw.vcf.processed.gz:GK010001
Using column name 'GK010002' for Indx19.fillmd.snp.raw.vcf.processed.gz:GK010002
Using column name 'GK010003' for Indx20.fillmd.snp.raw.vcf.processed.gz:GK010003
Using column name 'GK010004a' for Indx21.fillmd.snp.raw.vcf.processed.gz:GK010004a
Using column name 'GK010006' for Indx22.fillmd.snp.raw.vcf.processed.gz:GK010006
Using column name 'GK012201a' for Indx23.fillmd.snp.raw.vcf.processed.gz:GK012201a
Using column name 'GK012202' for Indx24.fillmd.snp.raw.vcf.processed.gz:GK012202
Using column name 'GK012203' for Indx25.fillmd.snp.raw.vcf.processed.gz:GK012203
Using column name 'GK012204' for Indx26.fillmd.snp.raw.vcf.processed.gz:GK012204
Using column name 'GK012301a' for Indx27.fillmd.snp.raw.vcf.processed.gz:GK012301a
Using column name 'GK012302' for Indx28.fillmd.snp.raw.vcf.processed.gz:GK012302
Using column name 'GK012303' for Indx29.fillmd.snp.raw.vcf.processed.gz:GK012303
Using column name 'GK012304' for Indx30.fillmd.snp.raw.vcf.processed.gz:GK012304
Using column name 'GK012305' for Indx31.fillmd.snp.raw.vcf.processed.gz:GK012305
Using column name 'GK016701a' for Indx32.fillmd.snp.raw.vcf.processed.gz:GK016701a
Using column name 'GK016702' for Indx33.fillmd.snp.raw.vcf.processed.gz:GK016702
Using column name 'GK016703' for Indx34.fillmd.snp.raw.vcf.processed.gz:GK016703
Using column name 'GK016704' for Indx35.fillmd.snp.raw.vcf.processed.gz:GK016704
Using column name 'GK016706' for Indx36.fillmd.snp.raw.vcf.processed.gz:GK016706
Using column name 'GK017801a' for Indx37.fillmd.snp.raw.vcf.processed.gz:GK017801a
Using column name 'GK017802' for Indx38.fillmd.snp.raw.vcf.processed.gz:GK017802
Using column name 'GK017803' for Indx39.fillmd.snp.raw.vcf.processed.gz:GK017803
Using column name 'GK017804' for Indx40.fillmd.snp.raw.vcf.processed.gz:GK017804
Using column name 'GK018701' for Indx41.fillmd.snp.raw.vcf.processed.gz:GK018701
Using column name 'GK018702' for Indx42.fillmd.snp.raw.vcf.processed.gz:GK018702
Using column name 'GK018703' for Indx43.fillmd.snp.raw.vcf.processed.gz:GK018703
Using column name 'GK018704a' for Indx44.fillmd.snp.raw.vcf.processed.gz:GK018704a
Using column name 'GK021002' for Indx45.fillmd.snp.raw.vcf.processed.gz:GK021002
Using column name 'GK021001' for Indx46.fillmd.snp.raw.vcf.processed.gz:GK021001
Using column name 'GK021003a' for Indx47.fillmd.snp.raw.vcf.processed.gz:GK021003a
Using column name 'GK021004' for Indx48.fillmd.snp.raw.vcf.processed.gz:GK021004
Using column name 'GK024101' for Indx49.fillmd.snp.raw.vcf.processed.gz:GK024101
Using column name 'GK024102' for Indx50.fillmd.snp.raw.vcf.processed.gz:GK024102
Using column name 'GK024103a' for Indx51.fillmd.snp.raw.vcf.processed.gz:GK024103a
Using column name 'GK024104' for Indx52.fillmd.snp.raw.vcf.processed.gz:GK024104
Using column name 'GK024105' for Indx53.fillmd.snp.raw.vcf.processed.gz:GK024105
Using column name 'GK027101a' for Indx54.fillmd.snp.raw.vcf.processed.gz:GK027101a
Using column name 'GK027102' for Indx55.fillmd.snp.raw.vcf.processed.gz:GK027102
Using column name 'GK027103' for Indx56.fillmd.snp.raw.vcf.processed.gz:GK027103
Using column name 'GK027104' for Indx57.fillmd.snp.raw.vcf.processed.gz:GK027104
Using column name 'GK027105' for Indx58.fillmd.snp.raw.vcf.processed.gz:GK027105
Using column name 'GK027701' for Indx59.fillmd.snp.raw.vcf.processed.gz:GK027701
Using column name 'GK034404' for Indx60.fillmd.snp.raw.vcf.processed.gz:GK034404
Using column name 'GK027703a' for Indx61.fillmd.snp.raw.vcf.processed.gz:GK027703a
Using column name 'GK027704' for Indx62.fillmd.snp.raw.vcf.processed.gz:GK027704
Using column name 'GK027705' for Indx63.fillmd.snp.raw.vcf.processed.gz:GK027705
Using column name 'GK028101' for Indx64.fillmd.snp.raw.vcf.processed.gz:GK028101
Using column name 'GK028102' for Indx65.fillmd.snp.raw.vcf.processed.gz:GK028102
Using column name 'GK028103a' for Indx66.fillmd.snp.raw.vcf.processed.gz:GK028103a
Using column name 'GK028104' for Indx67.fillmd.snp.raw.vcf.processed.gz:GK028104
Using column name 'GK029001' for Indx68.fillmd.snp.raw.vcf.processed.gz:GK029001
Using column name 'GK029002' for Indx69.fillmd.snp.raw.vcf.processed.gz:GK029002
Using column name 'GK029003' for Indx70.fillmd.snp.raw.vcf.processed.gz:GK029003
Using column name 'GK029004a' for Indx71.fillmd.snp.raw.vcf.processed.gz:GK029004a
Using column name 'GK029501' for Indx72.fillmd.snp.raw.vcf.processed.gz:GK029501
Using column name 'GK029502' for Indx73.fillmd.snp.raw.vcf.processed.gz:GK029502
Using column name 'GK029503a' for Indx74.fillmd.snp.raw.vcf.processed.gz:GK029503a
Using column name 'GK029504' for Indx75.fillmd.snp.raw.vcf.processed.gz:GK029504
Using column name 'GK034401a' for Indx76.fillmd.snp.raw.vcf.processed.gz:GK034401a
Using column name 'GK034402' for Indx77.fillmd.snp.raw.vcf.processed.gz:GK034402
Using column name 'GK034403' for Indx78.fillmd.snp.raw.vcf.processed.gz:GK034403
Using column name 'GK005005a' for Indx79.fillmd.snp.raw.vcf.processed.gz:GK005005a
Using column name 'GK035901a' for Indx80.fillmd.snp.raw.vcf.processed.gz:GK035901a
Using column name 'GK035902' for Indx81.fillmd.snp.raw.vcf.processed.gz:GK035902
Using column name 'GK035903' for Indx82.fillmd.snp.raw.vcf.processed.gz:GK035903
Using column name 'GK035904' for Indx83.fillmd.snp.raw.vcf.processed.gz:GK035904
Using column name 'GK035906' for Indx84.fillmd.snp.raw.vcf.processed.gz:GK035906
Using column name 'GK036501' for Indx85.fillmd.snp.raw.vcf.processed.gz:GK036501
Using column name 'GK036502' for Indx86.fillmd.snp.raw.vcf.processed.gz:GK036502
Using column name 'GK036503a' for Indx87.fillmd.snp.raw.vcf.processed.gz:GK036503a
Using column name 'GK036504' for Indx88.fillmd.snp.raw.vcf.processed.gz:GK036504
Using column name 'GK039101' for Indx89.fillmd.snp.raw.vcf.processed.gz:GK039101
Using column name 'GK039102' for Indx90.fillmd.snp.raw.vcf.processed.gz:GK039102
Using column name 'GK039103a' for Indx91.fillmd.snp.raw.vcf.processed.gz:GK039103a
Using column name 'GK039104' for Indx92.fillmd.snp.raw.vcf.processed.gz:GK039104
Using column name 'GK039601' for Indx93.fillmd.snp.raw.vcf.processed.gz:GK039601
Using column name 'GK039602' for Indx94.fillmd.snp.raw.vcf.processed.gz:GK039602
Using column name 'GK039603a' for Indx95.fillmd.snp.raw.vcf.processed.gz:GK039603a
Using column name 'GK039604' for Indx96.fillmd.snp.raw.vcf.processed.gz:GK039604
  • Use vcftools to convert to TPED/TFAM:
 ~ddiep/softwares/vcftools_0.1.9/bin/vcftools --vcf Merged_VCF --out UCLA_SZ_hg19_bissnp --plink-tped --recode
  • Change the chromosome names to numbers:
cp UCLA_SZ_hg19_bissnp.tfam Fixed.UCLA_SZ_hg19_bissnp.tfam
awk '{if($1 == "chrX") $1=23; if($1 == "chrY") $1=24; if($1 == "chrM") $1=26; print $0}' $f | sed 's/chr//g' > Fixed.$f
  • Use plink to filter (optional):
 ~ddiep/softwares/plink-1.07-x86_64/plink --noweb --geno 0.25 --tfile Fixed.UCLA_SZ_hg19_bissnp --recode --transpose --out UCLA_SZ_hg19_bissnp_filtered
  • plink reported some rsSNPs which were reported with multiple physical positions, we are not sure about them, use plink to remove those SNPs:
grep Duplicate UCLA_SZ_hg19_bissnp_filtered.log | awk -F " " '{print $6}' > snps_to_exclude
~ddiep/softwares/plink-1.07-x86_64/plink --noweb --exclude snps_to_exclude --tfile UCLA_SZ_hg19_bissnp_filtered --recode --transpose --out UCLA_SZ_hg19_bissnp_filtered-rmdup
  • Final TPED/TFAM:
UCLA_SZ_hg19_bissnp_filtered-rmdup.tped (9248 SNPs)
UCLA_SZ_hg19_bissnp_filtered-rmdup.tfam