Decontamination

Here, the script is carrying out alignment of the demultiplexed (and filtered) reads

Workflow

flowchart LR
subgraph Alignment
	al0[Filtered demultiplexed reads]
	al1{DNA or \ncDNA reads?}
	al2[Align vs DNA ref]
	al3[Align vs cDNA]
	al4[Flagstat BAM]
	al5[(Unaligned \nreads)]
	al6[de novo assembly workflow]
end
subgraph Functions in DNA_processing.py 
	dp1[align: index \nand align \nreads vs ref]
	dp2[DNA_align: calls \nalign func. \nto align DNA reads]
	dp3[cDNA_align: make/align \ntransciptome vs \ncDNA reads]
end
al0 --> al1
al1 --DNA--> al2
al1--cDNA-->al3
al2-.->al4
al2 --> al5
al3 --> al5
al5 --> al6
dp2 -..-> al2
dp3 -..-> al3
dp1 --> dp2
linkStyle 0 stroke: #008000
linkStyle 1 stroke: #008000
linkStyle 2 stroke: #008000
linkStyle 3 stroke: #008000
linkStyle 4 stroke: #008000
linkStyle 5 stroke: #008000
linkStyle 6 stroke: #008000
linkStyle 7 stroke: #8b4513
linkStyle 8 stroke: #ff1493
linkStyle 9 stroke: #14fff6

Loading

Alignment

• This is the step where the host/background or desired contaminant information is aligned against the filtered reads and eventually removed.
• For this, NanoMetaPipe.py calls the DNA_align and cDNA_align functions that are in the DNA_processing.py script
• Both DNA_align and cDNA_align functions also call a secondary function align that carries out the actual alignment command
• DNA_align: takes in:
  • the filtered/demultiplexed reads
  • the isolate name
  • the path to the stats directory
  • a temporary output folder that will be deleted after alignment
  • an alignment folder output
  • the number of threads
  • the path to the reference genome
• cDNA_align: take in:
  • the filtered/demultiplexed reads
  • the isolate name
  • if the user chooses to generate transcriptome:
    • the path to the forward and reverse reads to generate a transcriptome with
    • the path to the adapters
  • the output directory
  • the path to the NanoMetaPipe package's script folder
  • the path to the reference gff/transcriptome
  • the number of threads
  • the max memory to use
• If the user chooses to generate a transcriptome, cDNA_align uses a secondary script: cDNA_Processing.py to carry out the transcriptome assembly
  • Command to call the cDNA_Processing script

   cDNA_Processing.py -s reads -o path/to/output/cDNA_Processing -r path/to/forward_reads path/to/backwards_reads -a path/to/adapters -hr path/to/host_dna_reference -hg path/to/host_dna_gff -p threads -m max_memory -ur path/to/filterd_demultiplexed.fastq.gz

• Contamination removal from the metagenomic reads is achieved using minimap2. Example command is:

# index the reference/transcriptome
minimap2 -x map-ont -d path/to/reference.fasta.mmi path/to/reference.fasta
# align the reads against the genome/transcriptome
minimap2 -ax map-ont path/to/reference.fasta.mmi path/to/demultiplexed/reads.fastq -t threads | samtools view -@ threads -b - | samtools sort -@ threads -o path/to/output/alignments/readsVsRef.bam -

Extract unmapped reads

# stats using flagstat
samtools flagstat --threads path/to/output/alignments/readsVsRef.bam > path/to/output/stats/readsVsRef_FlagstatMappedVsRef_stats.txt
# use samtools to extract unaligned reads to carry forward
samtools view --threads threads -f 4 -b path/to/output/alignments/readsVsRef.bam > path/to/output/alignments/readsVsRef_unmapped.bam
# use bedtools to make fastq from the unmapped BAM file
bedtools bamtofastq -i path/to/output/alignments/readsVsRef_unmapped.bam -fq path/to/output/alignments/readsVsRef.bam

• After extracting the unmapped reads, de-deduplication is carried out using the filter_fastq_file function that is in the Tools.py
• The extracted unmapped reads is also quality assessed using NanoStat
• Extracted unmapped reads are taken forward for metagenome assembly and taxonomic classification