Coverage_depth_estimation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5347d63c",
   "metadata": {},
   "source": [
    "### Alignment\n",
    "The first script is the approach taken from D. Portugal-Barron. Unfortunately, I see to have missed a step here, or my approach, in which I substituted HISAT for Bowtie2, simply does not work. The second script is based on the following approach: https://www.metagenomics.wiki/tools/samtools/breadth-of-coverage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41851dd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "#define parameters:\n",
    "SAMPLENAME=\n",
    "PATTERNR1=name_pattern_R1\n",
    "PATTERNR2=name_pattern_R2\n",
    "GENOMEPATH=/path/to/genome/file/\n",
    "BOWTIE2_INDEXES=/path/to/indices\n",
    "READSPATH=/path/to/reads\n",
    "TEMPDIR=/temp/storage/dir/$SAMPLENAME/\n",
    "OUTDIR=/dir/for/results/$SAMPLENAME/\n",
    "#check whether indicated paths are correct\n",
    "\n",
    "[ -d $READSPATH ] || { echo 'Invalid path to mapping files, exiting' && exit; }\n",
    "[ -d $INDEXPATH ] || { echo 'Invalid path to contigs, exiting' && exit; }\n",
    "[ -d $GENOMEPATH ] || { echo 'Invalid path to contigs, exiting' && exit; }\n",
    "\n",
    "mkdir $OUTDIR\n",
    "mkdir $TEMPDIR\n",
    "for f in <name1> <name2> <name3>\n",
    "do\n",
    "\n",
    "#align reads to reference\n",
    "bowtie2 -x $BOWTIE2_INDEXES/\"$f\"_Index -1 $READSPATH/$PATTERNR1 -2 $READSPATH/$PATTERNR2 -S $TEMPDIR/\"$f\"_to_\"$SAMPLENAME\".sam\n",
    "#convert SAM to BAM\n",
    "samtools view -bS $TEMPDIR/\"$f\"_to_\"$SAMPLENAME\".sam > $TEMPDIR/\"$f\"_to_\"$SAMPLENAME\".bam\n",
    "#sort BAM file\n",
    "samtools sort $TEMPDIR/\"$f\"_to_\"$SAMPLENAME\".bam -o $OUTDIR/\"$f\"_to_\"$SAMPLENAME\".sorted.bam \n",
    "done; \n",
    "\n",
    "#then use samtools depth to estimate depth\n",
    "samtools depth $OUTDIR/*_to_\"$SAMPLENAME\".sorted.bam -o $OUTDIR/\"$SAMPLENAME\"_depth.txt\n",
    "#then use samtools coverage to estimate coverage\n",
    "samtools coverage *_to_\"$SAMPLENAME\".sorted.bam -o $OUTDIR/\"SAMPLENAME\"_coverage.txt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44f41b25",
   "metadata": {},
   "outputs": [],
   "source": [
    "#define parameters:\n",
    "SAMPLENAME=\n",
    "INDEXPATH=/path/to/indices\n",
    "READSPATH=/path/to/reads\n",
    "TEMPDIR=/temp/storage/dir/$SAMPLENAME/\n",
    "OUTDIR=/dir/for/results/$SAMPLENAME/\n",
    "\n",
    "[ -d $READSPATH ] || { echo 'Invalid path to mapping files, exiting' && exit; }\n",
    "[ -d $INDEXPATH ] || { echo 'Invalid path to contigs, exiting' && exit; }\n",
    "[ -d $GENOMEPATH ] || { echo 'Invalid path to contigs, exiting' && exit; }\n",
    "\n",
    "mkdir $OUTDIR\n",
    "mkdir $TEMPDIR\n",
    "#this makes a single file of the reads, allowing for a correct estimation of coverage and depth\n",
    "cat \"$READSPATH\"/\"$SAMPLENAME\"_*.fq > \"$TEMPDIR\"/\"$SAMPLENAME\"_reads.fq\n",
    "\n",
    "for f in Brevminu Durus Mcav SymClado SymKawa SymMicro SymNata SymNec SymPilo\n",
    "do\n",
    "#this creates the file needed to store the data\n",
    "echo 'These are the values for the depth of coverage at various target depths. First row is depth at 1X, then 5, 10, 15, 20, and 25. Last line is the target genome length.' > $OUTDIR/\"$f\"_\"$SAMPLENAME\".txt\n",
    "bowtie2 -x \"$INDEXPATH\"/\"$f\"_Index --no-unal -U \"$TEMPDIR\"/\"$SAMPLENAME\"_reads.fq -S \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sam -p $NSLOTS\n",
    "samtools view -bS \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sam > \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".bam\n",
    "samtools sort \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".bam > \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam\n",
    "samtools index \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam \n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=1 '$0>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=5 '$4>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=10 '$9>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=15 '$14>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=20 '$19>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "samtools mpileup \"$TEMPDIR\"/\"$f\"_\"$SAMPLENAME\".sorted.bam | awk -v X=25 '$24>=X' | wc -l >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "bowtie2-inspect -s \"$INDEXPATH\"/\"$f\"_index | awk '{ FS = \"\\t\" } ; BEGIN{L=0}; {L=L+$3}; END{print L}' >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "echo 'Depth of coverage at each given depth can be calculated by dividing the number of bases by the number of bases in the target genome.' >> \"$OUTDIR\"/\"$f\"_\"$SAMPLENAME\".txt\n",
    "done"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}