Skip to content

Instantly share code, notes, and snippets.

@tkrahn
Last active September 22, 2019 22:04
Show Gist options
  • Save tkrahn/612d21fb1691cecb11177e98aafb57b5 to your computer and use it in GitHub Desktop.
Save tkrahn/612d21fb1691cecb11177e98aafb57b5 to your computer and use it in GitHub Desktop.
This is my (simplified) pipeline to map Nanopore FastQ reads to a reference genome (here hg38). Note that this is certainly not the best way to use Nanopore results. It's just a quick check how the results look.
#!/bin/bash
START=$(date +%s.%N)
clear
# setup parameters
YSEQID=${PWD##*/}
# YSEQID="1234" # (the above command simply gets the name of the last segment of the current working directory)
NUM_THREADS=$(getconf _NPROCESSORS_ONLN)
echo "We can use ${NUM_THREADS} threads."
REF="/genomes/0/refseq/hg38/hg38.fa"
READS="fastq/${YSEQID}_nanopore.fastq.gz"
BAMFILE="${YSEQID}_minimap2_hg38.bam"
BAMFILE_SORTED="${YSEQID}_minimap2_hg38_sorted.bam"
# Pipeline commands:
minimap2 -t $NUM_THREADS $REF -ax map-ont $READS | \
samtools view -@ $NUM_THREADS -b -t $REF -o $BAMFILE -
samtools sort -@ $NUM_THREADS -T /var/tmp/sorted -o $BAMFILE_SORTED $BAMFILE
samtools index $BAMFILE_SORTED
samtools idxstats $BAMFILE_SORTED > ${BAMFILE_SORTED}.idxstats.tsv
# Delete no longer needed large file(s)
rm -f $BAMFILE # keep $BAMFILE_SORTED
# Easy tview file
echo "#!/bin/bash" > tview_${YSEQID}.sh
echo "samtools tview ${BAMFILE_SORTED} ${REF}" >> tview_${YSEQID}.sh
chmod a+x tview_${YSEQID}.sh
END=$(date +%s.%N)
DIFF=$(echo "$END - $START" | bc)
echo ${DIFF}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment