Changeset 5842

Timestamp:

Feb 25, 2020, 11:47:57 AM (3 years ago)

Author:

Nicklas Nordborg

Message:

References #1218: Implement MIPs alignment

Added functions for collecting metrics.

File:

• other/pipeline/trunk/mips_functions.sh (modified) (4 diffs)

other/pipeline/trunk/mips_functions.sh

@@ -192 +192 @@
 }
 
+# Get min_insert and max_insert from the Amplicons BED file
+# The values are needed by novoalign and collecting alignment metrics
+# Global parameters:
+#   $AmpliconsBed: Path to amplicons BED file
+#   $min_insert: Updated with value from the BED file
+#   $max_insert: Updated with value from the BED file
+function get_min_max_insert {
+  min_insert=$(awk 'NR == 1 || $3 - $2 < min {min = $3 - $2}END{print min - 1}' "${AmpliconsBed}")
+  max_insert=$(awk 'NR == 1 || $3 - $2 > max {max = $3 - $2}END{print max + 1}' "${AmpliconsBed}")
+}
+
+
 # Align with novoalign 
 # Parameters:
@@ -203 +215 @@
 #   $NovoSort: Path to novosort executable
 #   $NumThreads: Number of threads to use
+#   $min_insert: Calculated by get_min_max_insert
+#   $max_insert: Calculated by get_min_max_insert
 #   $AlignedBams: For returning paths to aligned BAM files (used for input in the merge step)
 function novoalign {
@@ -219 +233 @@
     -i PE ${min_insert}-${max_insert} \
     --amplicons ${AmpliconsBed} \
-    1 out/${prefix}.novo_coverage.bed \
+    1 out/${n}.novo_coverage.bed \
     -f tmp/${n}.umi.bam \
     > tmp/${n}.nosort.bam \
-    3> out/${prefix}.novo.log
+    3> out/${n}.novo.log
     
   ./stderrwrap.sh ${NovoSort} \
@@ -287 +301 @@
 }
 
+# Get PCR metrics for the concordant.bam file
+# Global parameters:
+#   $AmpliconsBed: Path to amplicons BED file
+#   $GenomeDict: Path to genome .dict file
+#   $GenomeFasta: Path to the genome fasta file
+#   $PcrMetricsOptions: Options for Picard
+function collect_pcr_metrics {
+  
+  local amplicon_set_name=$(basename ${AmpliconsBed} ".bed")
+  awk -F "\t" -v OFS="\t" '{print $1,$7,$8,$4,$5,$6}' ${AmpliconsBed} > amplicon_inserts.bed
+
+  ./stdwrap.sh ./picard2 BedToIntervalList \
+    -SEQUENCE_DICTIONARY ${GenomeDict} \
+    -INPUT amplicon_inserts.bed \
+    -OUTPUT amplicon_inserts.intervals \
+    >> pcr_metrics.out
+  
+  ./stdwrap.sh ./picard2 CollectTargetedPcrMetrics \
+    -INPUT out/concordant.bam \
+    -CUSTOM_AMPLICON_SET_NAME ${amplicon_set_name} \
+    -AMPLICON_INTERVALS amplicon_inserts.intervals \
+    -TARGET_INTERVALS amplicon_inserts.intervals \
+    -R ${GenomeFasta} \
+    ${PcrMetricsOptions} \
+    -PER_TARGET_COVERAGE out/concordant.per_amplicon_cov.txt \
+    -PER_BASE_COVERAGE out/concordant.per_base_cov.txt.gz \
+    -OUTPUT out/concordant.pcr_metrics.txt \
+    >> pcr_metrics.out
+}
+
+# Get alignment metrics for the concordant.bam file
+# Global parameters:
+#   $GenomeFasta: Path to the genome fasta file
+#   $max_insert: Calculated by get_min_max_insert
+#   $adapterN: Multiple adapter sequences defined earlier in this file
+function collect_alignment_metrics {
+
+  ./stdwrap.sh ./picard2 CollectAlignmentSummaryMetrics \
+    -INPUT out/concordant.bam \
+    -R ${GenomeFasta} \
+    -METRIC_ACCUMULATION_LEVEL null -METRIC_ACCUMULATION_LEVEL LIBRARY -METRIC_ACCUMULATION_LEVEL READ_GROUP \
+    -MAX_INSERT_SIZE ${max_insert} \
+    -ADAPTER_SEQUENCE null -ADAPTER_SEQUENCE ${adapter1} -ADAPTER_SEQUENCE ${adapter2} \
+    -ADAPTER_SEQUENCE ${adapter1rc} -ADAPTER_SEQUENCE ${adapter2rc} \
+    -OUTPUT out/concordant.alignment_summary_metrics.txt \
+    >> alignment_metrics.out
+
+}