# Can we add a column with the number of reads for each taxon next to relative\_abundance? **URL:** \ **Category:** MetaPhlAn\ **Created:** [July 3, 2023, 11:59am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577 "2023-07-03T11:59:40Z")\ **Posts on this page:** 8\ **Page:** 1
**Author:** ![Ameneh](https://avatars.discourse-cdn.com/v4/letter/a/ce73a5/32.png) [@Ameneh](https://forum.biobakery.org/u/Ameneh)\ **Post date:** [July 3, 2023, 11:59am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/1 "2023-07-03T11:59:41Z")
Hello, I am using **Metaphlan4** for fecal metagenomics. Although it works well, I would like to add the **exact number of reads** detected **for each taxon** in a column next to the **relative abundance** column. This is easy to calculate manually or in R/python since we are given the total number of reads in the sample, however would be equally useful to have it automatically calculated in the metaphlan output. Is there a particular parameter to get this that I am not aware of? I have read that some people used something like _“-t rel\_ab\_w\_read\_stats”_ but had a discrepancy in the counts obtained. Can you please help me? Best, Ameneh. ---
**Author:** ![aitor.blancomiguez](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/aitor.blancomiguez/32/86_2.png) [@aitor.blancomiguez](https://forum.biobakery.org/u/aitor.blancomiguez)\ **Post date:** [July 28, 2023, 1:55pm UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/2 "2023-07-28T13:55:24Z")
Hi @Ameneh MetaPhlAn 4 profiling is based on the mapping against species-specific markers genes, that only represent a part of the species genome. Thus the number of reads assigned to each species is just an estimation based on the avg coverage of the markers and the median length of the species’ genomes (what is calculated using “-t rel\_ab\_w\_read\_stats”). The EXACT number of reads is not possible to be calculated using MetaPhlAn ---
**Author:** ![Ameneh](https://avatars.discourse-cdn.com/v4/letter/a/ce73a5/32.png) [@Ameneh](https://forum.biobakery.org/u/Ameneh)\ **Post date:** [August 28, 2023, 1:43pm UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/3 "2023-08-28T13:43:19Z")
Hi @aitor.blancomiguez Thanks for your response. When I run Metaphlan 4, bowtie2out provide an text file as output with the number of reads. for example you can see this output: :#52373574 reads processed #SampleID Metaphlan\_Analysis #clade\_name NCBI\_tax\_id relative\_abundance additional\_species k\_\_Bacteria 2 100.0 k\_\_Bacteria|p\_\_Bacteroidetes 2|976 68.29302 … which means you have number of reads based on sequencing data you entered. But if there is something else that i am not aware of, please let me know, ---
**Author:** ![aitor.blancomiguez](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/aitor.blancomiguez/32/86_2.png) [@aitor.blancomiguez](https://forum.biobakery.org/u/aitor.blancomiguez)\ **Post date:** [September 1, 2023, 12:11pm UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/4 "2023-09-01T12:11:52Z")
Hi @Ameneh Yes, metaphlan will report the number of reads that were processed as input, but what I tried to say is that known the **exact** number of reads assigned to each taxa is not possible to assess via MetaPhlAn as it does not map against full genomes but just species-specific marker genes. The calc of the reads assigned to each taxa is a rough estimation based on the coverage of the markers and the median genome length of each species ---
**Author:** ![jihen](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/jihen/32/3342_2.png) [@jihen](https://forum.biobakery.org/u/jihen)\ **Post date:** [September 19, 2024, 7:57am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/5 "2024-09-19T07:57:51Z")
Hi @aitor.blancomiguez, Sorry for reviving this topic, but I have some questions regarding read counts in MetaPhlAn, especially after going through the 4.1 tutorial. I analyzed the file [SRS014476-Supragingival\_plaque.fasta.gz](https://github.com/biobakery/MetaPhlAn/releases/download/4.1.0/SRS014476-Supragingival_plaque.fasta.gz) and found the following: - Total reads: **20,000** (I counted this with `zcat SRS014476-Supragingival_plaque.fasta.gz | grep -c ">"`) - MetaPhlAn output: - Reads processed: **19,048** - Estimated reads mapped to known clades: **2,298** (using the `-t rel_ab_w_read_stats` flag) ### My Questions: 1. **From 20,000 to 19,048** : - You mention that “metaphlan will report the number of reads that were processed as input”, then how does the count decrease from **20,000** to **19,048**? Could the **952** reads that were not processed be due to them not being bacteria reads or possibly failing quality checks? 1. **From 19,048 to 2,298** : - I understand that the decrease to **2,298** may be influenced by the `-stat_q` parameter. When I tested different values: - `-stat_q 0.1` resulted in **4,659** estimated reads. - `-stat_q 0.0` resulted in **13,477** estimated reads. - If I set `-stat_q` to **0.0** , I expected to see the number of estimated reads equal to **19,048** since that would include all reads. Why is it still less than that? I appreciate your insights, and I hope I’m approaching this correctly! Thanks in advance! ---
**Author:** ![Claudia\_Mengoni](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/claudia_mengoni/32/2871_2.png) [@Claudia\_Mengoni](https://forum.biobakery.org/u/Claudia_Mengoni)\ **Post date:** [September 24, 2024, 9:07am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/6 "2024-09-24T09:07:19Z")
Hi @jihen to answer your questions: 1. The number of processed reads as you guessed is due to quality checks on reads, specifically based on minimum read length (default value is 70) 2. The 2,298 is the estimated number of reads that are assigned to a taxa by MetaPhlAn. The stat\_q parameter (default 0.20) is a parameter that is used to exclude upper and lower quartiles of markers detected for each taxa, so it will directly influence the estimation of reads. This is used to improve the accuracy of the profiling, reducing this to zero is not advised as it will increase sensitivity at the expense of specificity. The missing number of reads can be due again to the fact that this is an estimation, so it is unlikely it will correspond to the exact number of reads or to the fact that some reads may belong to unclassified taxa not present in the database. ---
**Author:** ![jihen](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/jihen/32/3342_2.png) [@jihen](https://forum.biobakery.org/u/jihen)\ **Post date:** [September 26, 2024, 12:45am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/7 "2024-09-26T00:45:20Z")
Hi @Claudia_Mengoni, thanks so much for your response! It really helped clarify my thoughts. ---
**Author:** ![jihen](https://yyz1.discourse-cdn.com/flex027/user_avatar/forum.biobakery.org/jihen/32/3342_2.png) [@jihen](https://forum.biobakery.org/u/jihen)\ **Post date:** [July 25, 2025, 7:47am UTC](https://forum.biobakery.org/t/can-we-add-a-column-with-the-number-of-reads-for-each-taxon-next-to-relative-abundance/5577/8 "2025-07-25T07:47:54Z")
Hi @Claudia_Mengoni, Sorry for necromancing this post again. I’m following up with another question after running MetaPhlAn 4.1 and comparing the stats with the input FASTQ files. Here’s the setup: - I input two paired-end FASTQ files (host removed), each with **99,320 reads**. - I verified that **all reads are ≥70 bp** using both `seqkit` and `awk`, e.g.: - So based on your previous response that MetaPhlAn filters out reads shorter than 70 bp, I think this show there will be no read drop as all reads are above the threshold. Then, I ran MetaPhlAn with this command: ```bash metaphlan RC2_host_rm_1.fq.gz,RC2_host_rm_2.fq.gz \ --input_type fastq -t rel_ab_w_read_stats \ --index mpa_vJun23_CHOCOPhlAnSGB_202307 \ --bowtie2db /home/bt2db/ \ --bowtie2out RC2.bowtie2.bz2 \ --nproc 8 -s RC2_metaphlan.sam -o RC2_profile_abs.txt ``` And I got this in the `RC2_profile_abs.txt` header: ```auto #198640 reads processed #estimated_reads_mapped_to_known_clades:82739 ``` But here’s where I get confused: - There are 99,320 reads are in each file, and all are ≥70 bp, while MetaPhlAn report **198,640 reads processed** (I guess this is because it added R1 + R2, so far so good), - Then only **82,739** reads are **estimated** to come from known clades (which as you said, some reads may belong to unclassified taxa or not present in the database, thus the number dropped), - However, when I export a sam file and do a flagstat, I saw only **4,944** reads appear as **mapped** in the SAM file: ### ❓ My questions: 1. **What is the relationship between:** 2. **Can I use `samtools flagstat` on the MetaPhlAn `.sam` file to infer anything meaningful about read mapping rate, coverage, or sequencing depth?** 3. **Similarly, can I use `samtools coverage` and `samtools depth` on the MetaPhlAn `.sam` file to infer anything meaningful about read mapping rate, coverage, or sequencing depth?** I’m trying to understand what each number represents in the pipeline and whether any of them can inform about read alignment quality or some statistics about sequencing coverage. * * * Thanks again for the great tool and support!