qc for research

alignment/bamqc.sh

@@ -35,6 +35,11 @@ shift $(($OPTIND -1))
 #    usage
 #fi
 
+if [[ -z $version ]]
+then
+    version='NA'
+fi
+
 source /etc/profile.d/modules.sh
 module load samtools/gcc/1.10 fastqc/0.11.8
 samtools flagstat ${sbam} > ${pair_id}.flagstat.txt
@@ -55,7 +60,7 @@ fi
 tmpdir=`pwd`
 if [[ $nuctype == 'dna' ]]; then
     module load bedtools/2.29.2 picard/2.10.3
-    bedtools coverage -sorted -g  ${index_path}/genomefile.txt -a ${bed} -b ${sbam} -hist > ${pair_id}.covhist.txt
+    bedtools coverage -a ${bed} -b ${sbam} -hist > ${pair_id}.covhist.txt
     grep ^all ${pair_id}.covhist.txt >  ${pair_id}.genomecov.txt
     perl $baseDir/calculate_depthcov.pl ${pair_id}.covhist.txt
     if [[ -z $skiplc ]]
@@ -63,13 +68,18 @@ if [[ $nuctype == 'dna' ]]; then
 	samtools view -@ $NPROC -b -L ${bed} -o ${pair_id}.ontarget.bam ${sbam}
 	samtools index -@ $NPROC ${pair_id}.ontarget.bam
 	samtools flagstat  ${pair_id}.ontarget.bam > ${pair_id}.ontarget.flagstat.txt
-	samtools view  -@ $NPROC -b -q 1 ${sbam} | bedtools coverage -sorted -hist -g ${index_path}/genomefile.txt -b stdin -a ${bed} > ${pair_id}.mapqualcov.txt
+	samtools view  -@ $NPROC -b -q 1 ${sbam} | bedtools coverage -hist -b stdin -a ${bed} > ${pair_id}.mapqualcov.txt
 	java -Xmx64g -Djava.io.tmpdir=${tmpdir} -XX:ParallelGCThreads=$NPROC -jar $PICARD/picard.jar EstimateLibraryComplexity BARCODE_TAG=RG I=${sbam} OUTPUT=${pair_id}.libcomplex.txt TMP_DIR=${tmpdir}
 	#java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectAlignmentSummaryMetrics R=${index_path}/genome.fa I=${pair_id}.ontarget.bam OUTPUT=${pair_id}.alignmentsummarymetrics.txt TMP_DIR=${tmpdir}
 	#samtools view  -@ $NPROC ${sbam} | awk '{sum+=$5} END { print "Mean MAPQ =",sum/NR}' > ${pair_id}.meanmap.txt
     fi
     #java -Xmx64g -Djava.io.tmpdir=${tmpdir} -jar $PICARD/picard.jar CollectInsertSizeMetrics INPUT=${sbam} HISTOGRAM_FILE=${pair_id}.hist.ps REFERENCE_SEQUENCE=${index_path}/genome.fa OUTPUT=${pair_id}.hist.txt TMP_DIR=${tmpdir}
-    perl $baseDir/sequenceqc_dna.pl -e ${version} -r $index_path ${pair_id}.genomecov.txt
+    if [[ $index_path/reference_info.pl ]]
+    then
+	perl $baseDir/sequenceqc_dna.pl -e ${version} -r $index_path ${pair_id}.genomecov.txt
+    else
+	touch ${pair_id}.genomecov.txt
+    fi
 else
     perl $baseDir/sequenceqc_rna.pl -e ${version} -r $index_path ${pair_id}.flagstat.txt
 fi

alignment/sequenceqc_dna.pl

@@ -122,7 +122,9 @@ foreach $sfile (@statfiles) {
 		 "Alignment_Status\t".$hash{status},"Alignment_Date\t".$hash{date},
 		 "File_Owner\t".$hash{fileowner},"Workflow_Version\t".$gittag),"\n";
   close OUT;
-  system(qq{cat $opt{refdir}\/reference_info.txt >> $prefix\.sequence.stats.txt});
+  if (-e "$opt{refdir}\/reference_info.txt") {
+      system(qq{cat $opt{refdir}\/reference_info.txt >> $prefix\.sequence.stats.txt});
+  }
   ##### END separateFilesPerSample ######
 }