add virus detection

Author: Shriram Bhosle

README.md

@@ -1,5 +1,12 @@
-# my_scripts
+##merge2files
 
 Script to compare and merge two files based on user specified matched columns.
 for usage:
 merge2files.pl -h
+
+##runVirusdetection
+* Tool to detect virus or bacterial sequences in tumour samples using GOTTCHA
+* Tracey Allen K. Freitas, Po-E Li, Matthew B. Scholz and Patrick S. G. Chain (2015) Accurate read-based metagenome characterization using a hierarchical suite of unique signatures, Nucleic Acids Research (DOI: 10.1093/nar/gkv180)
+* Download latest signature database for specific organism from : ftp://ftp.lanl.gov/public/genome/gottcha/GOTTCHA_database_v20150825/
+* Also needs lookup database GOTTCHA_lookup.tar.gz
+* runVirusdetection.sh 

merge2files.pl

@@ -9,28 +9,39 @@
 
 my ($options) = option_builder();
 
+
+if(!-s $options->{'f1'} && !-s $options->{'f2'}) {
+	warn "Unable to find files\n";
+	exit;
+} 
+
 open(my $fh1, $options->{'f1'}) || warn "print $!";
 open(my $fh2, $options->{'f2'}) || warn "print $!";
 
+
+
+
 my ($f1,$d1,$s1) = fileparse($options->{'f1'},qr/\.[^.]*/);
 my ($f2,$d2,$s2) = fileparse($options->{'f2'},qr/\.[^.]*/);
 
 my $col_f1=$options->{'c1'};
 my $col_f2=$options->{'c2'};
-$options->{'hdr'}=0 unless($options->{'hdr'}); 
+$options->{'hdr1'}=0 unless($options->{'hdr1'}); 
+$options->{'hdr2'}=0 unless($options->{'hdr2'}); 
 chomp $col_f1;
 chomp $col_f2;
 
+
 # create file handle for filed to write [ common columns and different columns
-open(my $fh_common, '>', $f1.'_'.$f2.'common'.$s1); 
+open(my $fh_common, '>', $options->{'o'}.'/'.$f1.'_'.$f2.'common'.$s1); 
 
 # create hash of both the files....
 
-my ($file1_hash)=get_file_hash($fh1,$col_f1,$options->{'hdr'},$f1);
-my ($file2_hash)=get_file_hash($fh2,$col_f2,$options->{'hdr'},$f2);
+my ($file1_hash)=get_file_hash($fh1,$col_f1,$options->{'hdr1'},$f1);
+my ($file2_hash)=get_file_hash($fh2,$col_f2,$options->{'hdr2'},$f2);
 
 #combine files
-combine_hash($file1_hash,$file2_hash,$fh_common,$options->{'hdr'});
+combine_hash($file1_hash,$file2_hash,$fh_common,$options->{'hdr1'});
 
 # sub to combine file data...
 sub combine_hash {
@@ -45,6 +56,7 @@ sub combine_hash {
   # write common rows
 	foreach my $key(keys %$hash_f1) {
 		next if $key eq 'header';
+		next if $key eq 'pre_header_lines';
 		if ($hash_f2->{$key}) {
 		  print $fh_combined $hash_f1->{$key}."\t".$hash_f2->{$key}."\n";
 		  $count++;
@@ -60,17 +72,17 @@ sub get_file_hash {
 	my $count=0;	
 	while (my $line=<$fh>) {	
 		$count++;
-		next if $line=~/^\s+\n/;
+		next if $line=~/^\n/;
 		$line=~s/\n//g;
 		my ($col_key) = get_col_key($line,$columns) if $count >= $header_row;
 		if($count <= $header_row ) {
 			if($count==$header_row) {
 			 $file_hash->{'header'}=$line;
+				print "Columns to compare from file: : $col_key\n" if defined $col_key; 
 			}
 			else {
 				$file_hash->{'pre_header_lines'}.="$line\n";
 			}
-			print "Columns to compare from file: : $col_key\n" if defined $col_key; 
 			next;
 		} 
 		$file_hash->{$col_key}="$line";
@@ -113,7 +125,9 @@ sub option_builder {
 					'f2|file2=s'  => \$opts{'f2'},
 					'c1|col1=s'  => \$opts{'c1'},
 					'c2|col2=s'  => \$opts{'c2'},
-					'hdr|header=s'  => \$opts{'hdr'},
+					'o|output=s'  => \$opts{'o'},
+					'hdr1|header1=s'  => \$opts{'hdr1'},
+					'hdr2|header2=s'  => \$opts{'hdr2'},
 					'v|version'  => \$opts{'v'},
 	);
 
@@ -129,7 +143,7 @@ =head1 NAME
 
 =head1 SYNOPSIS
 
-merge2files.pl [-h] -f1 -f2 -c1 -c2 [-hdr]
+merge2files.pl [-h] -f1 -f2 -c1 -c2 [-hdr1 -hdr2]
 
   Required Options:
     --help          (-h)  This message and format of input file
@@ -138,7 +152,9 @@ =head1 SYNOPSIS
     --file2  (-f2) file2 
     --col1  (-c1) comma separated column numbers fom file1  
     --col2  (-c2) comma separated columns numbers fom file2 
-    --header  (-hdr)  header row number 
+    --header1  (-hdr1)  header row number first file 
+    --header2  (-hdr2)  header row number second file
+    --output  (-o) output folder 
 
 =cut
 

runVirusDetection.sh

@@ -0,0 +1,51 @@
+#!/bin/bash
+set -ue
+##################
+##Author: sb43   #
+##Date:9/02/2016 #
+##################
+in_bam=${1:-}
+out_dir=${2:-}
+num_threads=${3:-}
+
+usage ()
+{
+  echo 'Usage : runVirusDetection.sh <BAM_FILE> <OUTDIR> <number_of_cpu>'
+  exit
+}
+
+if [[ -z "$in_bam" || -z "$out_dir" || -z "$num_threads" ]]
+then
+  usage
+fi
+
+mkdir -p $out_dir/tmp
+
+echo "Feteching unmapped reads......"
+#Unmapped read whose mate is mapped
+samtools view -u  -f 4 -F264  $in_bam > $out_dir/tmp/tmp1.bam &
+#mapped read whose mate is unmapped
+samtools view -u  -f 8 -F260  $in_bam > $out_dir/tmp/tmp2.bam &
+#Both reads unmapped
+samtools view -u  -f 12 -F256  $in_bam > $out_dir/tmp/tmp3.bam &
+# wait till all background processes finish
+wait
+# merge reads 
+echo "Merging unmapped reads......"
+samtools merge -u - $out_dir/tmp/tmp[123].bam | samtools sort -n - $out_dir/tmp/unmapped
+#create fastq
+echo "Creating fastq input for GOTTCHA script......"
+bedtools bamtofastq -i $out_dir/tmp/unmapped.bam -fq $out_dir/tmp/unmapped.fastq 
+
+echo " Running GOTTCHA script....."
+
+/software/CGP/external-apps/GOTTCHA/bin/gottcha.pl \
+--threads $num_threads \
+--mode all \
+--minQ 10 \
+--outdir $out_dir \
+--input $out_dir/tmp/unmapped.fastq \
+--database /lustre/scratch112/sanger/cgppipe/canpipe/test/ref/human/38/gottcha_db/database/GOTTCHA_VIRUSES_c5900_k24_u30_xHUMAN3x.strain
+
+echo "Completed virus detection using GOTTCHA....."
+echo "Check result file unmapped.gottcha.tsv .... "