# HG changeset patch # User geert-vandeweyer # Date 1411810697 14400 # Node ID 845a87ad254a8393d4b92ac2757b20f7113ccbdb # Parent a7e7f9b8f5387442d0c42803a6396678be356a53 re-upload, accidentally removed some files diff -r a7e7f9b8f538 -r 845a87ad254a DC_Genotyper.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DC_Genotyper.pl Sat Sep 27 05:38:17 2014 -0400 @@ -0,0 +1,530 @@ +#!/usr/bin/perl +############################# +## DEEP COVERAGE GENOTYPER ## +############################# +# version : 0.0.1 +# Principle: +# 1. get allele counts on all positions in specified targets (bed) using igvtools. Only SNPs !! +# 2. remove known dbsnp positions (bcf file) +# 3. Get distribution of background noise (pcr/sequencing errors), by modelling allele fractions as normal distributions. +# 4. Based on these distributions, check each position for significant change from the reference allele (based on allele fraction) +# 5. For abberant positions, check each alternate allele to see if it passes the background signal. +# 6. Generate VCF file. + + +################## +## LOAD MODULES ## +################## +use threads; +use threads::shared; +use Thread::Queue; +use Getopt::Std; + +#################### +## get paramaters ## +#################### +# t: target file +# b: bam file +# R: 2bit version of reference fasta. +# p: number of threads. +# s: dbsnp file +# m: minimal coverage (defaults 400x) +# P: ploidy +# a: outfile for allele distributions +# v: vcf file output. +getopts('t:b:R:p:s:m:P:v:a:', \%opts) ; + +## variables +my $twobit :shared; +my $igvgenome :shared; +if (!defined($opts{'R'})) { + die("Reference Genomes not specified\n"); +} +my @refgenomes = split(",",$opts{'r'}); +if (!-e $refgenomes[0]) { + die("'$refgenomes[0]' is not a valid file path."); +} +else { + $twobit = $refgenomes[0]; +} +if (!-e $refgenomes[1]) { + die("'$refgenomes[1]' is not a valid file path."); +} +else { + $igvgenome = $refgenomes[1]; +} + + +my $mincov :shared; +$mincov = 320; +if (defined($opts{'m'})) { + $mincov = $opts{'m'}; +} + +my $ploidy :shared; +if (defined($opts{'P'}) && $opts{'P'} =~ m/^\d+$/) { + $ploidy = $opts{'P'}; +} +else { + die("Ploidy (-P) was not specified or not an integer\n"); +} + + +if (defined($opts{'v'})) { + $outfile = $opts{'v'}; +} +else { + die("No output vcf-file specified.\n"); +} +if (!defined($opts{'a'})) { + die("No output file specified for distribution details\n"); +} +## create working dir. +my $rand = int(rand(10000)); +while (-d "/tmp/DC_Genotyper_$rand") { + $rand = int(rand(10000)); +} +my $wd :shared; +$wd = "/tmp/DC_Genotyper_$rand"; +system("mkdir '$wd'"); + + +my $snpfile :shared; +my $hassnp :shared; +$hassnp = 'NoDbSNP'; +$snpfile = ''; +if (defined($opts{'s'})) { + $snpfile = $opts{'s'}; + if (!-e $snpfile) { + die("'$snpfile' is not a valid file path."); + } + + my $mime = `file $snpfile`; + if ($mime !~ m/compressed/) { + print "$snpfile is not in compressed format. compressing & indexing the file now.\n"; + #print "... this takes a while\n"; + system("bgzip -c $snpfile > $wd/dbSNP.vcf.bgz"); + system("cd $wd/ && tabix -p vcf dbSNP.vcf.bgz"); + $snpfile = "$wd/dbSNP.vcf.bgz"; + } + elsif (!-e "$snpfile.tbi") { + print "tabix index file is missing for '$snpfile'. creating now.\n"; + ## check if I can write it out for future use + $snpfile =~ m/(.*)([^\/]+)$/; + my $d = $1; + if (-w $d) { + open OUT, ">$d/lock"; + flock(OUT,2); + system("cd $d && tabix -p vcf $snpfile"); + close OUT; + system("rm $d/lock"); + } + else { + system("cp $snpfile /$wd/dbSNP.vcf.bgz"); + system("cd $wd/ && tabix -p vcf dbSNP.vcf.bgz"); + $snpfile = "$wd/dbSNP.vcf.bgz"; + } + } + $hassnp = 'WithDbSNP'; +} + + +## 1. Get FASTA and prepare output hashes: +my $targets_one = Thread::Queue->new(); +my $targets_two = Thread::Queue->new(); +my $targets_three = Thread::Queue->new(); +open IN, $opts{'t'} or die("Could not open $opts{'t'} file for reading"); +if (-d "$wd/Fasta/") { + system("rm $wd/Fasta/*"); +} +else { + system("mkdir $wd/Fasta"); +} +## create the threads. +for (my $i = 1; $i<= $opts{'p'}; $i++) { + ${'thr'.$i} = threads->create('FetchFasta'); +} + +## enqueue the targets. +my %thash; +while () { + chomp; + my ($chr,$start,$stop,$name,$score,$strand) = split(/\t/,$_); + $targets_one->enqueue($_); + $targets_two->enqueue($_); + $targets_three->enqueue($_); + $thash{$chr}{$start} = $stop; +} +close IN; + +## end the threads. +for (my $i = 1; $i <= $opts{'p'}; $i++) { + $targets_one->enqueue(undef); +} + +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i}->join(); +} + +## load dbSNP inside target regions into shared structure. +########################################################## +my %dbsnp :shared; +if ($snpfile ne '') { + my $bcf = `which bcftools`; + chomp($bcf); + if ($bcf ne '') { + my $command = "bcftools query -f '\%CHROM\\t\%POS\\t\%REF\\t\%ALT\\t\%ID\\n' -R '".$opts{'t'}."' '$snpfile' > $wd/dbsnp.txt"; + system("$command"); + open IN, "$wd/dbsnp.txt"; + while () { + chomp; + my @p = split(/\t/,$_); + $dbsnp{$p[0].'-'.$p[1]} = $p[2].'-'.$p[3].'-'.$p[4]; + } + close IN; + } + else { + print "WARNING: BCFtools is not in the path. Skipping snp handling.\n"; + $snpfile = ''; + system("touch $wd/dbsnp.txt"); + } +} +else { + system("touch $wd/dbsnp.txt"); +} + +## now process the bam file. +mkdir "$wd/WIGS/"; +my $bam :shared; +$bam = $opts{'b'}; +my $bai = $bam.".bai"; +if (!-e $bai) { + print "BAI ($bai) missing for $bam : creating\n"; + system("samtools index $bam"); +} + +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i} = threads->create('CountAlleles'); +} +## end the threads. +for (my $i = 1; $i <= $opts{'p'}; $i++) { + $targets_two->enqueue(undef); +} + +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i}->join(); +} + +## generate the distributions. +############################## +my $alleles = Thread::Queue->new(); +my %all = ('A' => 1,'C' => 2,'G' => 3, 'T' => 4); +foreach(keys(%all)) { + $alleles->enqueue($_); + my $a = $_; + foreach(keys(%all)) { + if ($_ eq $a) { + next; + } + $alleles->enqueue($a.'-'.$_); + } +} +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i} = threads->create('GetDistribution'); +} +## end the threads. +for (my $i = 1; $i <= $opts{'p'}; $i++) { + $alleles->enqueue(undef); +} + +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i}->join(); +} + +## group distributions into one file +#################################### +my %map =('A' => 2,'C' => 3,'G' => 4, 'T' => 5); +open OUT, ">".$opts{'a'}; +print OUT "allele\tavg\tsd\n"; +foreach(keys(%map)) { + my $r = $_; + my $f = "$wd/model.$r.$mincov"."x.$hassnp.txt"; + open IN, "$f"; + my $a = ; + chomp($a); + #$dists{$r}{'avg'} = $a; + my $s = ; + chomp($s); + #$dists{$r}{'sd'} = $s; + close IN; + print OUT "$r\t$a\t$s\n"; + foreach(keys(%map)) { + if ($_ eq $r) { + next; + } + my $f = "$wd/model.$r-$_.$mincov"."x.$hassnp.txt"; + open IN, "$f"; + my $a = ; + chomp($a); + my $s = ; + chomp($s); + close IN; + print OUT "$r-$_\t$a\t$s\n"; + } +} +close OUT; + +## CALL SNPs +############ +# create the R script. +open R, ">$wd/CallSNPs.R"; +print R "args <- commandArgs(trailingOnly = TRUE)\n"; +print R "counts <- read.table(file=args[1],header=FALSE, as.is=TRUE)\n"; +print R "ploidy <- as.integer(args[3])\n"; +print R "chr <- args[2]\n"; +print R "snps <- read.table(file=args[5],header=FALSE,as.is=TRUE)\n"; +print R "colnames(snps) <- c('chr','pos','ref','alt','id')\n"; +print R "colnames(counts) <- c('pos','ref','A','C','G','T','TotalDepth')\n"; +print R "dists <- read.table(file='$wd/allelic_distributions.txt',header=TRUE,as.is=TRUE)\n"; +print R 'rownames(dists) = dists$allele'."\n"; +print R 'dists <- dists[,-1]'."\n"; +print R "vcf <- c()\n"; +print R "lower <- c()\n"; +print R "higher <- c()\n"; +print R "for (i in 1:(ploidy)) {\n"; +print R " lower[length(lower)+1] <- (2*i-1)/(2*ploidy)\n"; +print R " higher[length(higher)+1] <- (2*i+1)/(2*ploidy)\n"; +print R "}\n"; +print R "for (i in 1:nrow(counts)) {\n"; +print R " if (counts[i,'TotalDepth'] == 0) next\n"; +print R " # significantly different from reference?\n"; +print R " z <- ((counts[i,counts[i,'ref']]/counts[i,'TotalDepth']) - dists[counts[i,'ref'],'avg']) / dists[counts[i,'ref'],'sd']\n"; +print R " if (abs(z) > 3) {\n"; +print R " # test all alterate alleles to see which one is significant.\n"; +print R " for (j in c('A','C','G','T')) {\n"; +print R " if (j == counts[i,'ref']) next\n"; +print R " z <- ((counts[i,j]/counts[i,'TotalDepth']) - dists[paste(counts[i,'ref'],'-',j,sep=''),'avg']) / dists[paste(counts[i,'ref'],'-',j,sep=''),'sd']\n"; +print R " if (abs(z) > 3){\n"; +print R " filter <- 'PASS'\n"; +print R " phred <- round(-10*log(pnorm(-abs(z))),digits=0)\n"; +print R " if (phred > 9999) phred <- 9999\n"; +print R " frac <- counts[i,j]/counts[i,'TotalDepth']\n"; +print R " for (k in 1:ploidy) {\n"; +print R " if (frac >= lower[k] && frac < higher[k]) {\n"; +print R " sample <- paste(paste(paste(rep(0,(ploidy-k)),sep='',collapse='/'),paste(rep(1,k),sep='',collapse='/'),sep='/',collapse=''),':',counts[i,counts[i,'ref']],',',counts[i,j],sep='',collapse='')\n"; +print R " af <- k/ploidy\n"; +print R " break\n"; +print R " }\n"; +print R " }\n"; +print R " if (frac < lower[1]) {\n"; +print R " sample <- paste(paste(paste(rep(0,(ploidy-1)),sep='',collapse='/'),paste(rep(1,1),sep='',collapse='/'),sep='/',collapse=''),':',counts[i,counts[i,'ref']],',',counts[i,j],sep='',collapse='')\n"; +print R " af <- 1/ploidy\n"; +print R " filter <- 'LowFraction'\n"; +print R " }\n"; +print R " if (counts[i,'TotalDepth'] < $mincov) {\n"; +print R " filter <- 'LowCoverage'\n"; +print R " }\n"; +print R " info <- paste('DP=',counts[i,'TotalDepth'],';AF=',round(af,digits=5),';AR=',round(frac,digits=5),sep='')\n"; +print R " snpids <- which(snps\$chr == chr & snps\$pos == counts[i,'pos'])\n"; +print R " id <- '.'\n"; +print R " if (length(snpids) > 0) id <- snps[snpids[1],'id']\n"; +print R " vcf[length(vcf)+1] <- paste(chr,counts[i,'pos'],id,counts[i,'ref'],j,phred,filter,info,'GT:AD',sample,sep='\\t',collapse='')\n"; +print R " }\n"; +print R " }\n"; +print R " }\n"; +print R "}\n"; +print R "if (length(vcf) > 0) {\n"; +print R " write(file=args[4],paste(vcf,sep='\\n'))\n"; +print R "}\n"; +close R; +system("mkdir $wd/VCF/"); +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i} = threads->create('CallSNPs'); +} +## end the threads. +for (my $i = 1; $i <= $opts{'p'}; $i++) { + $targets_three->enqueue(undef); +} + +for (my $i = 1; $i <= $opts{'p'}; $i++) { + ${'thr'.$i}->join(); +} + +## BUILD FINAL VCF +open OUT, ">$outfile"; +print OUT "##fileformat=VCFv4.1\n"; +print OUT "##source=High_Ploidy_Genotyper_v.0.1\n"; +print OUT "##genome_reference=$twobit\n"; +if ($snpfile ne '') { + print OUT "##SNP_file=$snpfile\n"; +} +foreach(keys(%thash)) { + print OUT "##contig=\n"; +} +print OUT "##INFO=\n"; +print OUT "##INFO=\n"; +print OUT "##INFO=\n"; +print OUT "##FILTER=\n"; +print OUT "##FILTER=\n"; +print OUT "##FORMAT=\n"; +print OUT "##FORMAT=\n"; +print OUT "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tSAMPLE\n"; +close OUT; +@i = ( 1 .. 22,'X','Y','M' ); +foreach(@i) { + my $chr = "chr$_"; + foreach(sort {$a <=> $b} keys(%{$thash{$chr}})) { + my $v = "$wd/VCF/$chr.$_-".$thash{$chr}{$_}.".vcf"; + if (-e $v) { + system("cat '$v' >> '$outfile'"); + } + } +} + +## clean up +system("rm -Rf '$wd'"); + +sub FetchFasta { + while(defined(my $line = $targets_one->dequeue())) { + my ($chr,$start,$stop,$name,$score,$strand) = split(/\t/,$line); + # 2bit is zero based, non-including => decrease start by one + $startposition = $start - 1; + my $command = "twoBitToFa -seq=$chr -start=$startposition -end=$stop -noMask $twobit $wd/Fasta/$chr-$start-$stop.fasta"; + system($command); + } +} + +sub CountAlleles { + # local version of hashes + my $snp = \%dbsnp; + my %counts; + $counts{'A'} = ''; + $counts{'C'} = ''; + $counts{'G'} = ''; + $counts{'T'} = ''; + my %map =('A' => 1,'C' => 2,'G' => 3, 'T' => 4); + my %options; + foreach(keys(%map)) { + my $r = $_; + foreach(keys(%map)) { + if ($_ eq $r) { + next; + } + $options{$r.'-'.$_} = ''; + } + } + while (defined(my $line = $targets_two->dequeue())) { + $out = ''; + my ($chr,$start,$stop,$name,$score,$strand) = split(/\t/,$line); + ## get reference alleles + my %ref_alleles; + open FASTA, "$wd/Fasta/$chr-$start-$stop.fasta"; + my $head = ; + my $seq = ''; + while () { + chomp; + $seq .= $_; + } + close FASTA; + # this generates a hash of the reference alleles once, instead of substr-calls in every bam, on every iteration. + for (my $pos = 0; $pos < length($seq); $pos++) { + $ref_alleles{($pos+$start)} = substr($seq,$pos,1); + } + ## get counts. + my $target = "$chr:$start-$stop"; + my $command = "igvtools count -w 1 --bases --query '$target' '$bam' '$wd/WIGS/$chr-$start-$stop.wig' '$igvgenome' > /dev/null 2>&1"; + system($command); + open WIG, "$wd/WIGS/$chr-$start-$stop.wig"; + my $h = ; + $h = ; + $h = ; + my $target_counts = ''; + while () { + chomp; + #my ($pos, $a, $c, $g, $t , $n) = split(/\t/,$_); + my @p = split(/\t/,$_); + my $s = $p[1] + $p[2] + $p[3] + $p[4]; + $target_counts .= "$p[0]\t$ref_alleles{$p[0]}\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$s\n"; + ## skip positions with coverage < minimal coverage, and positions in dbsnp if specified (if not specified, snp hash is empty). + if ($s > $mincov && !defined($snp->{$chr.'-'.$p[0]})) { + ## for model of 'non-reference' + my $frac = $p[$map{$ref_alleles{$p[0]}}] / $s; + $counts{$ref_alleles{$p[0]}} .= $frac.','; + $out .= "$target\t$p[0]\t$ref_alleles{$p[0]}\t$p[1]\t$p[2]\t$p[3]\t$p[4]\n"; + ## for each of the options background models + foreach(keys(%map)) { + if ($_ eq $ref_alleles{$p[0]}) { + next; + } + $options{$ref_alleles{$p[0]}.'-'.$_} .= ($p[$map{$_}] / $s) .','; + } + + } + } + close WIG; + open OUT, ">>$wd/allcounts.$mincov"."x.$hassnp.txt"; + flock(OUT, 2); + print OUT $out; + close OUT; + open OUT, ">$wd/WIGS/$chr.$start-$stop.txt"; + print OUT $target_counts; + close OUT; + + } + foreach(keys(%counts)) { + open OUT, ">>$wd/counts_$_.$mincov"."x.$hassnp.txt"; + flock(OUT,2); + print OUT $counts{$_}; + close OUT; + } + foreach(keys(%options)) { + open OUT, ">>$wd/counts_$_.$mincov"."x.$hassnp.txt"; + flock(OUT,2); + print OUT $options{$_}; + close OUT; + } +} + +sub GetDistribution { + while (defined(my $allele = $alleles->dequeue())) { + system("sed -i 's/.\$//' '$wd/counts_$allele.$mincov"."x.$hassnp.txt'"); + open OUT, ">$wd/GetDistribution.$allele.R"; + print OUT "sample <- '$bam'\n"; + print OUT "nt <- '$allele'\n"; + #print OUT "pdf(file='$wd/Distribution.$allele.$mincov"."x.$hassnp.pdf',paper='a4')\n"; + print OUT "data <- scan(file='$wd/counts_$allele.$mincov"."x.$hassnp.txt',sep=',')\n"; + print OUT "nr <- length(data)\n"; + print OUT "avg <- mean(data)\n"; + print OUT "sdd <- sd(data)\n"; + #print OUT "if (avg > 0.5) {\n"; + #print OUT " x <- seq(0.8,1,length=1000)\n"; + #print OUT " y <- dnorm(x,mean=avg,sd=sdd)\n"; + #print OUT " plot(x,y,main='Distribution in sample $bam for nt $allele',xlab='Allelic Ratio',type='l',lwd=1)\n"; + #print OUT " abline(v=(avg-3*sdd),col='red')\n"; + #print OUT " text(0.81,max(y-0.5),paste(c('avg: ',avg,'\\nsd: ',sdd,'\\nnrDataPoints:', nr,'\\n$hassnp\\nMin.Cov:',$mincov),sep=' ',collapse=''),adj=c(0,1))\n"; + #print OUT "} else {\n"; + #print OUT " x <- seq(0,0.3,length=1000)\n"; + #print OUT " y <- dnorm(x,mean=avg,sd=sdd)\n"; + #print OUT " plot(x,y,main='Distribution in sample $bam for nt $allele',xlab='Allelic Ratio',type='l',lwd=1)\n"; + #print OUT " abline(v=(avg+3*sdd),col='red')\n"; + #print OUT " text(0.2,max(y-0.5),paste(c('avg: ',avg,'\\nsd: ',sdd,'\\nnrDataPoints:', nr,'\\n$hassnp\\nMin.Cov:',$mincov),sep=' ',collapse=''),adj=c(0,1))\n"; + #print OUT "}\n"; + #print OUT "dev.off()\n"; + print OUT "write(c(avg,sdd),file='$wd/model.$allele.$mincov"."x.$hassnp.txt',ncolumns=1)\n"; + close OUT; + system("cd $wd && Rscript GetDistribution.$allele.R >/dev/null 2>&1"); + } +} + + +sub CallSNPs { + while (defined(my $line = $targets_three->dequeue())) { + # split. + my ($chr,$start,$stop,$name,$score,$strand) = split(/\t/,$line); + my $file = "$wd/WIGS/$chr.$start-$stop.txt"; + my $ofile = "$wd/VCF/$chr.$start-$stop.vcf"; + system("cd $wd && Rscript CallSNPs.R '$file' '$chr' '$ploidy' '$ofile' '$wd/dbsnp.txt'"); + } + +} + diff -r a7e7f9b8f538 -r 845a87ad254a DC_Genotyper.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/DC_Genotyper.xml Sat Sep 27 05:38:17 2014 -0400 @@ -0,0 +1,83 @@ + + + + R_3_0_2 + samtools + tabix + blat_server + perl_module_threads + perl_module_threads_shared + perl_module_Thread_Queue + igvtools + + DC_Genotyper.pl + -t "$targets" + -b "$bamfile" + -R "${ref.fields.path}" + -p "\${GALAXY_SLOTS:-4}" + #if $dbsnp.source == "history": + -s "${dbsnp.ownFile}" + #else + -s "${dbsnp.indices.fields.path}" + #end if + -m $mincov + -P $ploidy + + -a $output1 + -v $output2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +**What it does** + + 1. get allele counts on all positions in specified targets (bed) using igvtools. Only SNPs !! + 2. remove known dbsnp positions (bcf file) + 3. Get distribution of background noise (pcr/sequencing errors), by modelling allele fractions as normal distributions. + 4. Based on these distributions, check each position for significant change from the reference allele (based on allele fraction) + 5. For abberant positions, check each alternate allele to see if it passes the background signal. + 6. Generate VCF file. + + +**Information** + +This tools is created by Geert Vandeweyer. It is a very early version with several limitations. Current limitations are : no support for indels, no plotting of the noise-models, incorrect syntax in for multi-allelic sites in the VCF file. + +Any feedback is welcome. + + + diff -r a7e7f9b8f538 -r 845a87ad254a readme.rst --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/readme.rst Sat Sep 27 05:38:17 2014 -0400 @@ -0,0 +1,18 @@ +BACKGROUND: + +DC_Genotyper stands for Deep-Coverage Genotyper, and is aimed at detecting low fraction SNPs (no indels) in high-ploidy (or pooled) samples with very high coverage. It is being developed at the University of Antwerp by Geert Vandeweyer. + +METHOD: + +DC_Genotyper generates a background noise distributions on a per-sample basis, and uses these distrubutions to detect non-reference sites. For non-reference sites, Allele-specific distributions are used to estimate if an allele surpasses the background signal (e.g. A_to_G has different distribution as A_to_C, also reflected in Tr/Tv ratios). + +LIMITATION: + +This is a very early version with several limitations. Current limitations are : no support for indels, no plotting of the noise-models, incorrect syntax in for multi-allelic sites in the VCF file. + +Any feedback is welcome + + +INSTALLATION: + +After installation, complete the dbsnp.loc and dc_genotyper_indexes.loc files. The DC_genotyper_indexes.loc file has a specific format, specifying two index files per genome in a single line. Make sure you use comma to seperate these. DCG supports multithreading, but keep tests have shown that using more than 6-8 threads will lead to I/O bottlenecks. diff -r a7e7f9b8f538 -r 845a87ad254a tool_dependencies.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tool_dependencies.xml Sat Sep 27 05:38:17 2014 -0400 @@ -0,0 +1,28 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + +