Skip to content

Classification confirmation

Jump to bottom
Damilola Oresegun edited this page Dec 19, 2022 · 7 revisions

This is another optional secondary pipeline that is to be carried out after the completion of the NanoMetaPipe main pipeline. As previously stated, the main pipeline concludes in the taxonomic classification by kraken2 and subsequent generation of a visual interactive chart. To confirm the classification determined by Kraken2, a secondary confirmatory workflow using a hidden markov model (HMM) profile was developed. Importantly, this secondary workflow was constructed and tested using small subunit rRNA however, it should fundamentally work for other sequences that can be targeted using a HMM approach. Other similar sequence targets could include things like the merozoite surface protein (MSP) or other species/genus specific targets.

Workflow

flowchart TD
subgraph PreChecks
direction LR
a{User: run or create?}
a1[run]
a2[create]
a2a{candidates file provided?}
a2a2[end]
a1a{hmm profile provided?}
a1b[end]
a1a1{profile indexed?}
a1a1a[delete indexes]
a3[continue]
a4{zipped inputs?}
a4a[unzip]
a5{FASTQ input?}
a5a[convert to FASTA]
a6[start analysis]
end
a --> a1
a --> a2
a1 --> a1a
a1a --yes--> a1a1
a1a --no--> a1b
a1a1 --yes--> a1a1a
a1a1 --no--> a3
a2 --> a2a
a2a --yes--> a3
a2a --no--> a2a2
a1a1a --> a3
a3 --> a4
a4 --yes--> a4a
a4a --> a6
a4 --no-->a5
a5 --yes--> a5a
a5a --> a6
a5 --no--> a6
linkStyle 0 stroke: #adff2f
linkStyle 1 stroke: #adff2f
linkStyle 2 stroke: #adff2f
linkStyle 3 stroke: #adff2f
linkStyle 4 stroke: red
linkStyle 5 stroke: #adff2f
linkStyle 6 stroke: #adff2f
linkStyle 7 stroke: #adff2f
linkStyle 8 stroke: #adff2f
linkStyle 9 stroke: red
linkStyle 10 stroke: #adff2f
linkStyle 11 stroke: #adff2f
linkStyle 12 stroke: #adff2f
linkStyle 13 stroke: #adff2f
linkStyle 14 stroke: #adff2f
linkStyle 15 stroke: #adff2f
linkStyle 16 stroke: #adff2f
linkStyle 17 stroke: #adff2f
subgraph Run Mode
direction LR
r1[align reads to HMM]
r1a[index HMM profile]
r1b[align reads to indexed HMM profile]
r1c[Output]
r1d[[RawMatches.txt]]
r1e[[PerSequenceMatches.txt]]
r1f[[Short_PerHit.txt]]
r2[extracted aligned reads]
r2a[extract aligned readIDs]
r2b[extract full read sequence]
rbc[[IsolateID_SSU_reads.fasta]]
r3[search against SILVA database]
r3a{reads de-replicated?}
r3a1[de-replicate reads]
r3b{database de-replicated?}
r3b1[de-replicate database]
r4[search reads against database]
r5[[VSEARCH_clusters.tab]]
r6[[VSEARCH_BLAST_alignedReads.txt]]
a6 --> r1
r1 --> r1a
r1a --> r1b
r1b --> r1c
r1c -.-> r1d
r1c -.-> r1e
r1c -.-> r1f
r1e --> r2
r2 --> r2a
r2a --> r2b
r2b -.-> rbc
rbc --> r3
r3 --> r3a
r3a --yes--> r3b
r3a --no--> r3a1
r3a1 --> r3b
r3b --yes--> r4
r3b --no--> r3b1
r3b1 --> r4
r4 -.-> r5
r4 -.-> r6
end
linkStyle 18 stroke: #008000
linkStyle 19 stroke: #008000
linkStyle 20 stroke: #008000
linkStyle 21 stroke: #008000
linkStyle 22 stroke: #996600
linkStyle 23 stroke: #996600
linkStyle 24 stroke: #996600
linkStyle 25 stroke: #008000
linkStyle 26 stroke: #008000
linkStyle 27 stroke: #008000
linkStyle 28 stroke: #008000
linkStyle 29 stroke: #008000
linkStyle 30 stroke: #008000
linkStyle 31 stroke: #008000
linkStyle 32 stroke: #008000
linkStyle 33 stroke: #008000
linkStyle 34 stroke: #008000
linkStyle 35 stroke: #008000
linkStyle 36 stroke: #008000
linkStyle 37 stroke: #996600
linkStyle 38 stroke: #996600
subgraph Create Mode
direction LR
c1[Download Candidate Sequences]
c1a[[DownloadedCandidates.fasta]]
c1b{File empty?}
c1b1[Alert user. End]
c1b2[Continue]
c2[Multiple sequence alignment]
c2a[Align candidate sequences with eachother]
c2a1[[my_prefix_candidates_aligned.aln]]
c2b[remove poorly aligned regions]
c2b1[[my_prefix_candidates_aligned_clean.fasta]]
c3[Construct HMM profile]
c3a[[my_prefix.hmm]]
c4[Pass HMM profile to the 'run mode' workflow]
a6 --> c1
c1 -.-> c1a
c1a --> c1b
c1b --yes--> c1b1
c1b --no--> c1b2
c1b2 --> c2
c2 --> c2a
c2a -.-> c2a1
c2a1 --> c2b
c2b -.-> c2b1
c2b1 --> c3
c3 -.-> c3a
c3a --> c4
c4 --> r1
end
linkStyle 39 stroke: #008000
linkStyle 40 stroke: #996600
linkStyle 41 stroke: #008000
linkStyle 42 stroke: #008000
linkStyle 43 stroke: #008000
linkStyle 44 stroke: #008000
linkStyle 45 stroke: #008000
linkStyle 46 stroke: #996600
linkStyle 47 stroke: #008000
linkStyle 48 stroke: #996600
linkStyle 49 stroke: #008000
linkStyle 50 stroke: #996600
linkStyle 51 stroke: #008000
linkStyle 52 stroke: #008000
Loading

Arguments

Required arguments

  • -f or --fasta: Full path to the your query fasta file. These are reads to be mapped against the provided/generated HMM profile. Example: path/to/Unaligned_reads.fasta
  • -m or --mode: Choose which mode to run the script. In create mode, the script will create a HMM profile but will require a candidates file (-c) to do this. In run mode, the script will require a previously made HMM profile (-p) to run the rest of the script. Choices include: create or run. Example: -m create
  • -o or --output: Full path to the directory to save generated outputs. Example: path/to/project/TaxonomicClassification
  • -db or --database: Full path to the directory that holds the SILVA database

Optional Arguments

  • -c or --candidates: Full path to a .txt file of a list of accession IDs for candidate sequences to generate the HMM profile. This argument should only be used if you want to generate a HMM profile. Example: path/to/candidates_Accession.txt
  • -p or --profile: Full path to a previously generated HMM profile made from DNA sequences. This should be used only in run mode.
  • -pp or --profile_prefix: Use only in create mode. Use this option to provide a specific prefix name to save the created HMM profile. Default is MyHMMProfile
  • -s or --isolate: Give a name for the isolate. If not chosen, the script will attempt to find the isolate name from either filename or the directory name. Default is Isolate
  • -incE or --inclusion_evalue: State the inclusion expectation value threshold. E-values above the value given are excluded from the outputs. Default is 0.01 which is an expectation of 1 false positive in 100 query searches
  • -E or --e_value: State the expectation value for false positives. This option works to allow some noise through and for the use to determine the most appropriate hiti based on their data needs. Default is 5 i.e. no more than 5 false positives may be returned per query sequence
  • -t or --threads: Give a number of threads to use. Default is 2
  • -vt or --vsearch_threshold: State a number between 0 and 1 to act as a similarity threshold for vsearch global search. Note: the higher the threshold value, the more stringent and thus the less the number of reads returned in the BLAST output

To run

To run the script with a previously generated profile:

# activate the correct conda environment name 'hmmENV'
conda activate hmmENV
# run the SSU script with a previously generated profile
python ./SSU_rRNA.py -f full/path/to/reads.fasta -m run -o path/to/output -db path/to/SILVA/database

To first generate a HMM profile before confirmation is carried out:

# activate the correct conda environment name 'hmmENV'
conda activate hmmENV
# run the SSU script to create a new HMM profile
python ./SSU_rRNA.py -f full/path/to/reads.fasta -m create -o path/to/output -db path/to/SILVA/database

Explanation

This workflow is essentially either creating or using a HMM profile to confirm the identification of a target of interest that was determined in the kraken2 taxonomic classification. In the case of the development of this pipeline, the target of interest was the SSU rRNA of Plasmodium. As such, the SILVA database downloaded only contained the genus of Plasmodium. Firstly the script checks the input arguments to determine if the user has chosen the create or run modes. Where the user chooses to create, checks are done to see if the user provides the candidates.txt file before proceeding. If the user chooses to run the script, checks are done to see if the hmm profile is provided and if it has been previoiusly indexed. Additionally, checks are done to determine if the inputs are zipped and subsequently unzips if this is the case. Furthermore, where the input is a FASTQ file, it is automatically converted to FASTA.

Run mode

  • inputted reads are aligned against the provided HMM profile using the hmmAln function. The function
    • indexes the HMM profile
    • aligns the reads against the HMM profile using the nhmmscan function of the hmmer package
    • outputs aligned reads information in three files:
      • RawMatches.txt
      • PerSequenceMatches.txt
      • Short_PerHit.txt
  • After this, reads identified to align to the HMM profile are extracted using the getReads function. The function
    • extracts the readIDs of reads that were identified from nhmmscan
      • saves identified readIDs to isolate_List_of_Hit_Reads.txt
    • uses the extracted readIDs to extract the corresponding full length reads using the fastaExtract function
      • extracted reads are saved in isolate_SSU_reads.fasta
  • The extracted reads are then mapped against the SILVA data using the vsrch function. The function
    • checks if the reads have been de-replicated. If not, then de-replication is done
    • checks if the database has been de-replicated. If not then, de-replication is done
    • outputs two files:
      • VSEARCH_clusters_UC.tab
      • VSEARCH_BLAST_alignedReads.txt

Create mode

  • The downloadCandidates function is used to carry out downloaded provided candidate sequences from NCBI. The function
    • uses the secondary script Download_NCBI_Accession_to_Fasta.py to download the full length sequences corresponding to the candidate accession codes
    • saves the downloaded candidates to a FASTA file: DownloadedCandidates.fasta
    • checks if the DownloadedCandidates.fasta is not empty and alerts the user if it is empty
  • carries out a multiple sequence alignment (MSA) of the DownloadedCandidates.fasta using the alnMuscle function
    • the resulting alignment is saved as my_prefix_candidates_aligned.aln
    • trimal is used to remove poorly aligned and spurious sequences
    • cleaned alignment file is saved as my_prefix_candidates_aligned_clean.fasta
  • constructs the HMM profile using hmmbuild with the provided prefix
  • passes the create HMM profile to the run mode workflow described above

Known Issue

As stated, this workflow was initially developed to investigate the SSU_rRNAs for Plasmodium however, this was not entirely successful. While the pipeline was able to identify reads which align to the SSU database derived from SILVA, it was found that many of the reads identified were not unique to Plasmodium spp.. This does not seem to be an error in the pipeline but rather in the target sequence. It is our thought that the SSU_rRNA target was not unique enough to be appropriate for the dataset and a target sequence that is unique to just Plasmodium spp. would be considerable more appropriate. As such, this workflow is still viable however, this is down to the user to determine if it is appropriate for their use-case, research question and dataset.

Clone this wiki locally