Human and phiX read-removal with kraken2

--- tags: cov-irt --- # Human and phiX read-removal with kraken2 This took about 20 minutes to set up as run below. [toc] ## Creating conda environment ```bash conda create -n kraken2 -c conda-forge -c bioconda -c defaults kraken2=2.0.9beta conda activate kraken2 ``` ## Setting up kraken2 database with human and phiX genomes --- > **NOTE** > This was updated on 11-Sept-2020, built the same way except for one change, adding the `--no-masking` flag to the `kraken2-build --download library` and `kraken2-build --add-to-library` commands. The reason for this is discussed at the end of this page [here](https://hackmd.io/Qd0XKEsUR1ebOnOPubqbrg#Why-building-with-no-masking). Roughly following along with [here](https://github.com/DerrickWood/kraken2/wiki/Manual#custom-databases). Downloading NCBI taxonomy info needed (takes like 5 minutes): ```bash mkdir kraken2_human_and_phiX_db kraken2-build --download-taxonomy --db kraken2_human_and_phiX_db/ ``` Downloading human reference (takes ~1 minute as run here): ```bash kraken2-build --download-library human --db kraken2_human_and_phiX_db/ \ --threads 40 --no-masking ``` Downloading and adding phiX genome to this library: ```bash curl -L https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/819/615/GCF_000819615.1_ViralProj14015/GCF_000819615.1_ViralProj14015_genomic.fna.gz | gunzip > phiX.fa kraken2-build --add-to-library phiX.fa --db kraken2_human_and_phiX_db/ --no-masking # removing genome fasta rm phiX.fa ``` Building database (takes ~7 minutes as run here): ```bash kraken2-build --build --db kraken2_human_and_phiX_db/ --threads 40 ``` Removing intermediate files: ```bash kraken2-build --clean --db kraken2_human_and_phiX_db/ ``` ### Download database as built on 11-Sept-2020 (LATEST) 3GB compressed, ~4.3GB uncompressed. Can be downloaded with the following: ```bash curl -L -o kraken2_human_and_phiX_db.tar.gz https://ndownloader.figshare.com/files/24658262 tar -xzvf kraken2_human_and_phiX_db.tar.gz ``` ### Download database as built on 19-June-2020 > This is the older one that was built with masking, and I don't think we should use anymore. See [below](https://hackmd.io/Qd0XKEsUR1ebOnOPubqbrg#Why-building-with-no-masking) for why the change. 2.8GB compressed, ~4GB uncompressed. Can be downloaded with the following: ```bash #### NOTE: the date is on the tar.gz file, but will not be part of the directory name once it is unpacked curl -L -o kraken2_human_and_phiX_db_19_June_2020.tar.gz https://ndownloader.figshare.com/files/23276654 tar -xzvf kraken2_human_and_phiX_db_19_June_2020.tar.gz ``` ## Example command filtering out human and phiX reads Getting tiny example data: ```bash curl -L -o sample-1.fq.gz https://ndownloader.figshare.com/files/23237460 ``` Performing filtering: ```bash kraken2 --db kraken2_human_and_phiX_db/ --threads 20 \ --output sample-1-kraken2-out.txt --report sample-1-kraken2-report.txt \ --unclassified-out sample-1-filtered-reads.fq sample-1.fq.gz # then compressing the output if wanted gzip sample-1-filtered-reads.fq ``` ## Tracking info on reads removed per sample Here's one way we can make a table that has some info like percent reads removed from each sample. As written, it assumes output files named like above ({sample-ID}-kraken2-out.txt), and would need a single-column text file with the unique sample IDs as input. For an example, here's copying the output from the above to just make a second output for example purposes: ```bash cp sample-1-kraken2-out.txt sample-2-kraken2-out.txt cp sample-1-kraken2-report.txt sample-2-kraken2-report.txt ``` Making a sample list input file: ```bash printf "sample-1\nsample-2\n" > sample-list.txt ``` Pasting and running this next codeblock will generate the bash script to do the summary: ```bash cat << 'EOF' > summarizing-kraken2-read-removal.sh #!/usr/bin/env bash # making sure we aren't adding to one already started rm -rf building.tmp for SAMPLE_ID in $(cat $1); do total_fragments=$(wc -l ${SAMPLE_ID}-kraken2-out.txt | cut -f 1 -d " ") fragments_retained=$(grep -w -m 1 "unclassified" ${SAMPLE_ID}-kraken2-report.txt | cut -f 2) perc_removed=$(printf "%.2f\n" $(echo "scale=4; 100 - ${fragments_retained} / ${total_fragments} * 100" | bc -l) ) printf "${SAMPLE_ID}\t${total_fragments}\t${fragments_retained}\t${perc_removed}\n" >> building.tmp done # combining cat <( printf "Sample_ID\tTotal_reads_before\tTotal_reads_after\tPercent_reads_removed\n" ) building.tmp > kraken2-read-removal-summary.tsv && rm building.tmp EOF ``` Making script executable: ```bash chmod +x summarizing-kraken2-read-removal.sh ``` Now that script can be put wherever we want, but the sample input list needs to point to where the kraken2 output files are stored. So probably easiest to run it in the directory that has the kraken2 output files, then that list can just be the file names, as in this example: ```bash ./summarizing-kraken2-read-removal.sh sample-list.txt ``` Which gives us this summary table: ``` Sample_ID Total_reads_before Total_reads_after Percent_reads_removed sample-1 2 1 50.00 sample-2 2 1 50.00 ``` --- # Why build with no masking? First, "masking" is the process of masking low-complexity regions so they aren't seen/considered by whatever process we are going to put them through. There is a commonly used program called [dustmasker](https://www.ncbi.nlm.nih.gov/IEB/ToolBox/CPP_DOC/lxr/source/src/app/dustmasker/) that does this that's used by NCBI, kraken2, centrifuge, and others. While working on some read-based classifier databases, running one of our samples that already had been run against the original 19-June-2020 human and phiX kraken2 db to remove human reads, I was getting a ton of human reads classified. I was working specifically with [this file](https://osf.io/hg7p4/) if curious, and about 1 million out of 2 million reads were getting classified as human – again, this is *after* going through the human-read removal process. This was happening with both centrifuge and a different kraken2 database I had access to (that both held human, bacteria, archaea, and fungi), and when I pulled out the reads getting classified as human and ran them through BLAST, they sure enough were coming back clearly human. So, after driving myself crazy for quite a long, wonderful time trying to figure out WTF was going on, it all came down to this masking process. kraken2 does masking by default, and generally it's a good idea. If we want to *classify* things, I'd use masking to build that database. In that case we don't want ambiguous regions giving us false positives. But if we want to *identify and remove something* like is the case here, I don't think masking is helpful. Aside from these being inherently low-information sequences by definition, there are only two possibilities for reads that would match to an un-masked, low-complexity fragment of the human genome: 1. They only match human, in which case we would want them removed anyway 2. They are ambiguous enough that they came from something else but they still match at least the human genome too, and potentially other things, in which case we wouldn't want to trust their classification and would want them removed anyway In both of those cases, we don't need or want them. So for building a database like this that's sole purpose is to remove things that match the human or phiX genome, I think it's better to build the reference DB without masking.

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.