parsing virus status, msisensor-pro

alignment/combine_idxstats.pl → alignment/filter_viral_idxstats.pl

@@ -8,11 +8,16 @@ use Getopt::Long;
 
 #my $cwd = getcwd();
 
+my $pair_id = "";
+GetOptions ('pairid|p=s' => \$pair_id);
+if($pair_id eq ""){$pair_id = "all";}
+
 my @allFiles = @ARGV;
 chomp(@allFiles);
+
 my %data;
 my @samples;
-open OUT, ">viral.info" or die $!;
+open OUT, ">$pair_id\.viral.seqstats.txt" or die $!;
 
 my %virus;
 $virus{'NC_001538.1'}="BK polyomavirus";

alignment/viralalign.sh → alignment/virusalign.sh

@@ -9,12 +9,13 @@ usage(){
 }
 
 OPTIND=1 #Reset OPTIN
-while getopts :b:r:p:h opt
+while getopts :b:r:p:fh opt
 do
   case $opt in
     b) bam=$OPTARG;;
     r) ref=$OPTARG;;
     p) pairid=$OPTARG;;
+    f) filter=1;;
     h) usage;;
   esac
 done
@@ -35,6 +36,7 @@ fi
 reffa=${ref}/idt_virus_reference.fa
 source /etc/profile.d/modules.sh
 module load bwa/intel/0.7.17 picard/2.10.3 samtools/1.6 
+baseDir="`dirname \"$0\"`"
 
 samtools view -@ 8 -b -u -F 2 ${bam} |samtools sort -n - >unmapped.bam
 java -Djava.io.tmpdir=./ -Xmx4g -jar $PICARD/picard.jar SamToFastq I=unmapped.bam FASTQ=unmapped.R1.fastq SECOND_END_FASTQ=unmapped.R2.fastq UNPAIRED_FASTQ=unmapped.unpaired.fastq
@@ -45,4 +47,7 @@ samtools view -h -F 256 -b out.sam -o out.bam
 samtools sort out.bam -o ${pairid}.viral.bam
 samtools index ${pairid}.viral.bam
 samtools idxstats ${pairid}.viral.bam >${pairid}.viral.idxstats.txt
-samtools flagstat ${pairid}.viral.bam >${pairid}.viral.flagstat.txt
+if [[ $filter == 1 ]]
+then
+    perl $baseDir/filter_viral_idxstats.txt -p ${pairid} ${pairid}.viral.idxstats.txt
+fi

variants/msisensor.sh

@@ -9,13 +9,14 @@ usage() {
   exit 1
 }
 OPTIND=1 # Reset OPTIND
-while getopts :r:l:n:b:p:h opt
+while getopts :r:l:n:c:b:p:h opt
 do
     case $opt in
         r) index_path=$OPTARG;;
         p) pair_id=$OPTARG;;
         b) sbam=$OPTARG;;
         n) normal=$OPTARG;;
+	c) capture=$OPTARG;;
         h) usage;;
     esac
 done
@@ -32,11 +33,17 @@ fi
 source /etc/profile.d/modules.sh	
 export PATH=/project/shared/bicf_workflow_ref/seqprg/bin:$PATH
 
+bedopt=''
+if [[ -n $capture ]]
+then
+    bedopt="-e $capture"
+fi
+   
 if [[ -n $normal ]]
 then
-    msisensor2 msi -d ${index_path}/microsatellites.list -n $normal -t $sbam -o ${pair_id}.msi
     
+    msisensor-pro msi -d ${index_path}/microsatellites.list -n $normal -t $sbam -o ${pair_id}.msi $bedopt
 else
     # -M ${index_path}/msi_tumoronly_model 
-    msisensor2 msi -d ${index_path}/microsatellites.list -t $sbam -o ${pair_id}.msi
+    msisensor2 msi -d ${index_path}/microsatellites.list -t $sbam -o ${pair_id}.msi $bedopt
 fi