Dinh/Dinh 2012/NOTES/2012-9-22
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