annotate DC_Genotyper.pl @ 13:08c4fb95c3a9 draft

changed dep to samtools .19 to include bcftools

#!/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: reference genome files for twobit and IGV.
# 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 (<IN>) {
	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 (<IN>) {
			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 = <IN>;
	chomp($a);
	#$dists{$r}{'avg'} = $a;
	my $s = <IN>;
	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 = <IN>;
		chomp($a);
		my $s = <IN>;
		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=<ID=$_,assembly=hg19,species=\"Homo Sapiens\">\n";
}
print OUT "##INFO=<ID=DP,Number=1,Type=Integer,Description=\"Total Depth\">\n";
print OUT "##INFO=<ID=AF,Number=1,Type=Float,Description=\"Allele Frequency\">\n";
print OUT "##INFO=<ID=AR,Number=1,Type=Float,Description=\"Allelic Ratio\">\n";
print OUT "##FILTER=<ID=LowFraction,Description=\"Allelic Fraction under 1/2*$ploidy\">\n";
print OUT "##FILTER=<ID=LowCoverage,Description=\"Total Depth is lower than threshold of $mincov\">\n";
print OUT "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">\n";
print OUT "##FORMAT=<ID=AD,Number=2,type=Integer,Description,\"Allelic Depth\">\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 = <FASTA>;
		my $seq = '';
		while (<FASTA>) { 
			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 = <WIG>;
		$h = <WIG>;
		$h = <WIG>;
		my $target_counts = '';
		while (<WIG>) {	
			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'");
	}

}