Hi all,
I recently wrote a python script that can normalize ASVs with 16S rRNA gene copy number. Hope it works well and suggestions are welcome!
Introduction:
Normalize sequences by 16S rRNA gene copy number (GCN) based on rrnDB database (version 5.6). The script matches the taxa of sequences with the rrnDB-5.6_pantaxa_stats_NCBI.tsv
file, starting from the lowest rank. If a match is found, the mean of GCN for the taxon is assigned; if not, the script will try to match a higher rank until the highest rank is met. All the unassigned sequences are assumed to have one GCN.
Note that the mean column in the rrnDB-5.6_pantaxa_stats_NCBI.tsv
is, according to the rrnDB manual, calculated from the means of the pan-taxa of immediate lower rank. Therefore, the mean of GCN might be different from the rrndb online search result. For example, the “mean” of GCN for bacteria is 2.02 in the downloading tsv file, whereas the mean of GCN for all the bacterial taxa is 5.0 if you search rrnDB online database.
Setting path:
First clone the repository using command below:
git clone https://github.com/Jiung-Wen/16S_copy_num_normalize.git
To set up the path, open .bash_profile
(macOS) or .bashrc
(Linux) with any text editor you prefer:
vim ~/.bash_profile
In .bash_profile
, append the following command:
export PATH="/YOUR_PATH/16S_copy_num_normalize/copy_num_normalize/:$PATH"
Save and close the file. Restart your terminal or using command below to apply the change immediately :
source ~/.bash_profile
Usage:
We assume that you have installed and activated QIIME2 environment.
copy_num_normalize.py --table table.qza --taxonomy taxonomy.qza -d silva -o output_file_name
-
--table
PATH - path of QIIME2 artifactFeatureTable[Frequency]
-
--taxonomy
PATH - path of QIIME2 artifactFeatureData[Taxonomy]
-
-d
STRING - database used for sequence annotation {silva, greengenes} -
-o
PATH - path of output directory and file name (If path is not gave, output files will save to current directory.)
Running example:
We use artifacts from QIIME2’s “Moving Pictures” tutorial as test files. Use the following commands to download the files.
# DADA2 output artifact:
wget https://docs.qiime2.org/2019.10/data/tutorials/moving-pictures/table-dada2.qza
# Taxonomic analysis output artifact:
wget https://docs.qiime2.org/2019.10/data/tutorials/moving-pictures/taxonomy.qza
We can normalize the FeatureTable using the command below:
copy_num_normalize.py --table table-dada2.qza \
--taxonomy taxonomy.qza \
-d greengenes \
-o table-dada2
The outputs would be a GCN normalized artifact table-dada2_copy_number_normalized.qza
of type FeatureTable[Frequency]
and a .txt file table-dada2_16S_rRNA_copy_number.txt
that indicates the GCN for each sequence.
Now you can perform analyses as you usually do in QIIME2 with the GCN-normalized FeatureTable.