source: other/pipeline/trunk/mips_functions.sh @ 5838

# $Id: mips_functions.sh 5838 2020-02-24 14:45:46Z nicklas $
# Helper functions for MIPS alignment. This script can't be used by itself.
# It has been designed to be included in a main script that also need to
#  set several global variables. For details see each function definition.
#
# Copy FASTQ and RG files from/to given folders
# Parameters: 
#   $1: Folder with FASTQ and RG files
#   $2: Folder to copy the files to          
function copy_fastq_and_rg {
  local srcFolder=$1
  local destFolder=$2
  if [ ! -d "${srcFolder}" ]; then
    echo "Can't find FASTQ data folder ${srcFolder}" 1>&2
    exit 1
  fi
  cp ${srcFolder}/*.fastq.gz $2
  cp ${srcFolder}/*.rg $2
}

# Adapter sequences
# Trimmomatic adapter sequences:
# "For 'Palindrome' clipping, a matched pair of adapter sequences must be provided.
#  The sequence names should both start with 'Prefix', and end in '/1' for the
#  forward adapter and '/2' for the reverse adapter. The part of the name between
# 'Prefix' and '/1' or '/2' must match exactly within each pair."
# Checking the adapter fasta files provided by Trimmomatic, it seems as the
# Prefix/1 adapter is the i5 adapter read from the flowcell surface (i.e. NOT
# as it will appear in read2 when reading through a short fragment)
# Prefix/2 adapter is the i7 adapter, also as read from the flowcell surface 
# MIPS BDC i5 adapter sequence after the index for trimmomatic: CATACGAGATCCGTAATCGGGAAGCTGAAG
# MIPS BDC i7 adapter sequence after the index for trimmomatic: ACACGCACGATCCGACGGTAGTGT
adapter1=CATACGAGATCCGTAATCGGGAAGCTGAAG
adapter2=ACACGCACGATCCGACGGTAGTGT
# Picard CollectAlignmentSummaryMetrics adapter sequences:
# I haven't figured out in what direction it wants the adapters,
# for simplicity, I use both directions
# Reverse complement of MIPS BDC i5 adapter: CTTCAGCTTCCCGATTACGGATCTCGTATG
# Reverse complement of MIPS BDC i7 adapter: ACACTACCGTCGGATCGTGCGTGT
adapter1rc=CTTCAGCTTCCCGATTACGGATCTCGTATG
adapter2rc=ACACTACCGTCGGATCGTGCGTGT
# Novoalign adapter sequences:
# Adapter1 is the i7 adapter sequence as it appears in read1 when reading
# through a short fragment
# Adapter2 is the i5 adapter sequence as it appears in read2 when reading
# through as short fragment
# MIPS BDC read1 (i7) adapter (before i7 index): ACACTACCGTCGGATCGTGCGTGT
# MIPS BDC read2 (i5) adapter (before i5 index): CTTCAGCTTCCCGATTACGGATCTCGTATG
# Use 12 first bases of each for novoalign (if using -a)
novo_adapter1=ACACTACCGTCG
novo_adapter2=CTTCAGCTTCCC

# Create an adapter file that can be used with Trimmomatic for adapter trimming
# Parameters:
#   $1: Name of the file to create
function create_adapter_fa {
  local faFile=$1
  echo -e ">PrefixPE/1\n${adapter1}\n>PrefixPE/2\n${adapter2}" > ${faFile}
}

# Run two Trimmomatic steps. First trim adapters and then quality
# Parameters:
#   $1: Prefix to FASTQ files "_R?.fastq.gz" is added automatically
#   $2: Prefix for generated files (counter)
# Global parameters:
#   $JAVA: Path to Java exectable
#   $TrimmomaticJAR: Path to Trimmomatic JAR file
#   $NumThreads: Number of threads to use
#   $TrimOptionsAdapter: Options for adapter trimming
#   $TrimOptionsQual: Options for quality trimming
#   $nAfterTrim: Counter for number of FASTQ files with data remaining after trimming
function run_trimmomatic {
  local prefix=$1
  local n=$2
  
  ## Adapter trimming
  ./stdwrap.sh ${JAVA} -jar ${TrimmomaticJAR} PE \
    -threads ${NumThreads} \
    -phred33 \
    -trimlog tmp/${n}.adaptrim.log \
    fastq/${prefix}_R1.fastq.gz fastq/${prefix}_R2.fastq.gz \
    tmp/${n}.adaptrim.1.fastq tmp/${n}.adaptrim.u.1.fastq tmp/${n}.adaptrim.2.fastq tmp/${n}.adaptrim.u.2.fastq \
    ${TrimOptionsAdapter} \
    >> trimmomatic.out
    
  if [ ! -s "tmp/${n}.adaptrim.1.fastq" ]; then
    rm -f tmp/${n}.adaptrim.*.fastq
    echo "${prefix}: No paired end reads surviving adapter trimming"
    return
  fi
    
  ## Quality trimming
  ./stdwrap.sh ${JAVA} -jar ${TrimmomaticJAR} PE \
    -threads ${NumThreads} \
    -phred33 \
    -trimlog tmp/${n}.qualtrim.log \
    tmp/${n}.adaptrim.1.fastq tmp/${n}.adaptrim.2.fastq \
    tmp/${n}.qualtrim.1.fastq tmp/${n}.qualtrim.u.1.fastq tmp/${n}.qualtrim.2.fastq tmp/${n}.qualtrim.u.2.fastq \
    ${TrimOptionsQual} \
    >> trimmomatic.out
    
  if [ ! -s "tmp/${n}.qualtrim.1.fastq" ]; then
    rm -f tmp/${n}.qualtrim.*.fastq
    echo "${prefix}: No paired end reads surviving quality trimming"
    return
  fi
  
  nAfterTrim=$((nAfterTrim + 1))
}

# Generate statistics from Trimmomatic log files
# Parameters:
#   $1: Input file name pattern (*.log files from Trimmomatic)
#   $2: Output file name
function trimlog_stats {
  local input=$1
  local output=$2
  
  echo -e "Trim\tCount" > ${output}
  if ls $input > /dev/null 2>&1; then
    awk '{print $ (NF)}' ${input} | sort -n | uniq -c | awk -v OFS="\t" '{print $2,$1}' >> ${output}
  fi
}

# Read RG file and convert to Picard2 syntax.
# Parameters:
#   $1: Path to RG file
# Returns:
#   Converted RG value
function read_and_fix_rg {
  local rg=$1
  
  read -r DATA < ${rg}
  
  DATA=${DATA/RG=/-READ_GROUP_NAME }
  DATA=${DATA/SM=/-SAMPLE_NAME }
  DATA=${DATA/LB=/-LIBRARY_NAME }
  DATA=${DATA/PU=/-PLATFORM_UNIT }
  DATA=${DATA/DT=/-RUN_DATE }
  DATA=${DATA/PL=/-PLATFORM  }
  DATA=${DATA/CN=/-SEQUENCING_CENTER }
  
  echo ${DATA}
}

# Run Picard to convert FASTQ files to BAM file
# Parameters:
#   $1: Prefix to RG file ".rg" is added automatically
#   $2: Prefix for generated files (counter)
function fastq_to_bam {
  local prefix=$1
  local n=$2
  
  if [ ! -f tmp/${n}.qualtrim.1.fastq ]; then
    return
  fi
  
  local rg=$(read_and_fix_rg "fastq/${prefix}.rg")
  
  ./stdwrap.sh ./picard2 FastqToSam \
    -SORT_ORDER queryname \
    -QUALITY_FORMAT Standard \
    ${rg} \
    -FASTQ tmp/${n}.qualtrim.1.fastq \
    -FASTQ2 tmp/${n}.qualtrim.2.fastq \
    -OUTPUT tmp/${n}.u.bam \
    >> fastqtosam.out
}

# Run FgBio to annotate the BAM file with UMI information
# Parameters:
#   $1: Prefix to UMI file "_UMI.fastq.gz" is added automatically
#   $2: Prefix for generated files (counter)
# Global parameters:
#   $JAVA: Path to Java exectable
#   $FgBioJAR: Path to FgBio JAR file
function annotate_umi {
  local prefix=$1
  local n=$2
  
  if [ ! -f tmp/${n}.u.bam ]; then
    return
  fi
  
  ./stdwrap.sh ${JAVA} -jar ${FgBioJAR} AnnotateBamWithUmis \
    --fail-fast true -t RX \
    -i tmp/${n}.u.bam \
    -f fastq/${prefix}_UMI.fastq.gz \
    -o tmp/${n}.umi.bam \
    >> annotatebam.out

}

# Align with novoalign 
# Parameters:
#   $1: Prefix of sample (used in output files)
#   $2: Prefix for generated files (counter)
# Global parameters:
#   $NovoAlign: Path to novoalign exectable
#   $NovoAlignOptions: Options for novoalign
#   $NovoIndex: Path to index file 
#   $AmpliconsBed: Path to amplicons BED file
#   $NovoSort: Path to novosort executable
#   $NumThreads: Number of threads to use
#   $AlignedBams: For returning paths to aligned BAM files (used for input in the merge step)
function novoalign {
  local prefix=$1
  local n=$2

  if [ ! -f tmp/${n}.umi.bam ]; then
    return
  fi

  ./stderrwrap.sh ${NovoAlign} \
    -c ${NumThreads} \
    ${NovoAlignOptions} \
    -d ${NovoIndex} \
    --tags MC,ZP \
    -i PE ${min_insert}-${max_insert} \
    --amplicons ${AmpliconsBed} \
    1 out/${prefix}.novo_coverage.bed \
    -f tmp/${n}.umi.bam \
    > tmp/${n}.nosort.bam \
    3> out/${prefix}.novo.log
    
  ./stderrwrap.sh ${NovoSort} \
    -c ${NumThreads} -t . -i \
    -o tmp/${n}.novo.bam \
    tmp/${n}.nosort.bam \
    3> tmp/${n}.sort.log

  AlignedBams="${AlignedBams} -INPUT tmp/${n}.novo.bam"
}

# Merge BAM files
# Global parameters:
#   $AlignedBams: BAM files to merge (created by 'novoalign')
function merge_bam {

  ./stdwrap.sh ./picard2 MergeSamFiles \
    -SORT_ORDER coordinate -ASSUME_SORTED true \
    ${AlignedBams} \
    -OUTPUT tmp/novo.bam \
    > mergebam.out
}