qichao1984 / ncyc Goto Github PK

As mentioned in the title, I want to annotate the N genes in some MAGs, if I can use the ORF predicted from prodigal? And I want to know the difference between the two databases. Before use, are they must be unzipped?

Hello

I am interested in how many comammox amoA genes are in my samples, but I don`t know whether this database contain these genes. If contains, at this time can I separate them from other amoA_A genes ?

thanks
hu

Hello! Can I assign more threads to process? Its using only one processor.

Can output from different runs be combined or is the abundance relative to all other present samples?

I have 3 large assemblies and would rather be allowed to run the in parallel. Is this allowed or must they be run together?

-d : specify the directory where your fasta/fastq (or gzipped) files are located.

Hi,
Thanks for the great database. I used NCYC on my MAG's contig and got the results. one thing I want to make sure is to get the particular gene sequence used for annotation. For say,
G01_000000067461 0fd91084f91d2de2731de1f1356dd7d6 nirB

In the above result, it says G01_000000067461 corresponds to nirB gene, but in reality, G01_000000067461 is a part of a contig/scaffold representing multiple genes. So how do I can get the particular gene sequence from G01_000000067461 used by NCYC for nirB.

Similarly, the diamond output file does not have any headers, and it's a bit confusing to what some of the number and characters correspond to.

Diamond output
G01_000000067461 0fd91084f91d2de2731de1f1356dd7d6 80.6 813 155 2 69690 67258 1 812 0.0e+00 1329.3

Hi All,

I'm afraid blastall -d flag only receives custom-built or NCBI standard database, e.g. nr, refseq.
Raw fasta files must be built into a custom blast database in prior for blastp.
But in the line 77 of NCycProfiler.PL:
system("blastall -p blastp -d data/NCyc_100.faa -i $file -o $out -m 8 -e 1e-4 -b 1 -v 1")
attempts to directly assign the faa file to the -d flag. Has anyone made this successful?
BTW, the use of blastall is also discontinued from legacy blast suite. For users using newer versions of blast, modify the script to use blastp/blastn instead.

In addition, the error message of the above error is confusing and not documented. Since no return value checked after the system call, the fail of blastall is not caught. The script keeps execution and ended up in an error message from cp instead: no such file or directory of some_sample.blast (due to no output from blastall).

-m diamond has no this issue. Not checked with usearch mode.

Hello and thank you for this very useful tool for finding N-cycling related genes!

I have a question regarding the output of using the NCyc_95 vs NCyc_100(_2019Jul) databases: how is it that the 95 database can produce more hits for a gene than the 100 database? To provide some background, I am running the program on a set of MAGs and wanted to see how the output would differ using the two different databases. I have multiple MAGs where there are a greater number of hits for genes resulting from using the 95 database compared to the 100 database. Why does this happen?

I am curious about this because it would seem that any of the representative sequences belonging to the 95 database should also be present in the 100 database? Am I misunderstanding how these two databases are created?

I am running the tool using Diamond and am happy to share any code or resulting output of the program as you may need! Thank you for your time and advice!

Hi,
Is it possible to provide the accession IDs of the sequences used in NCyc_100.faa.gz?
Thanks,
Ruben

Hi,
Thanks for this great work! I'm interested in using it.
Could you please provide a short example on how to use the NCycProfiler.PL script for profiling metagenomic sequence data (e.g., fasta of gene sequences).
Thanks,
Ruben

I assume you parse the diamond output to somehow generate the count tables? I am wondering how you do this as I would like to know on what contigs the genes were found so I can generate a pseudo abundance of that gene based on mapping depth to that contig.

I wish to do this as a count is only telling me how contigs with that gene on exist but hides the potential biological significant of how abundant that contig may be (i.e. how abundant that bacteria it came from it)

Hi,

I have a similar problem to issue #11. I am running NCyc on two samples as a test. These files are called sample31merged.fastq.gz (127189806 reads)and sample32merged.fastq.gz (191219006 reads).
I run the following command:

perl NCycProfiler.PL -d ~/Dorset_mesocosms/ProjectLIMS19618/CLEAN_READS2/MERGED_READS/ncyctest/ -m diamond -f fastq.gz -s nucl -si ~/Dorset_mesocosms/ProjectLIMS19618/CLEAN_READS2/MERGED_READS/ncyctest/samples.tsv -o ~/nitro_genes.tsv

After running I get the following:

Deallocating buffers... [0.043s]
Deallocating queries... [0.004s]
Loading query sequences... [0s]
Closing the input file... [0.018s]
Closing the output file... [0.132s]
Closing the database file... [0.01s]
Deallocating taxonomy... [0s]
Total time = 5829.35s
Reported 3850528 pairwise alignments, 3850528 HSPs.
3850528 queries aligned.
The host system is detected to have 1621 GB of RAM. It is recommended to increase the block size for better performance using these parameters : -b12 -c1
was not found in /home/bgregs/Dorset_mesocosms/ProjectLIMS19618/CLEAN_READS2/MERGED_READS/ncyctest/samples.tsv, please check!

Diamond appears to successfully run generating sample31merged.diamond and sample32merged.diamond files but then I get the above error and there is no further output. This seems to be an issue with reading the sample information file.
I've tried a few different formats with the sample info file. I have tried it with (samples.tsv) or without (samples2.tsv) the fastq.gz extension as recommended in issue #14 and without the full file path (sampleinfo.tsv) but I get the same error regardless. I have attached these files as .txt files as I can't upload tsv to GitHub:

sampleinfo.txt
samples2.txt
samples.txt

The only changes I have made to the script is path to diamond:

##Please specify where your prefered database searching tool locates
##NCBI blast ftp://ftp.ncbi.nlm.nih.gov/blast/executables/legacy.NOTSUPPORTED/2.2.26/blast-2.2.26-x64-linux.tar.gz
my $blast = "/bin/blast/bin/blastall";
my $formatdb = "/bin/blast/bin/formatdb";
##diamond https://github.com/bbuchfink/diamond/releases
my $diamond = "/usr/local/bin/diamond2";
##usearch https://www.drive5.com/usearch/download.html
my $usearch = "~/bin/usearch8.1.1861_i86linux32";

Any advice to sort out this issue would be appreciated.

Hi!

I am trying to make a blast db with the NCyc100.faa.gz file (I also tried with the 95 and got the same error) so that I can blast it outside of the perl executable you provide. I basically unzip the file and then run makeblastdb, However I can't make the database because i get the following error:

`BLAST Database creation error: Error: Duplicate seq_ids are found:
LCL|39946.BGIOSGA039077-PA

`
Has this issue happened before?
Thank you!

Hello,

I made functional annotation associated with nitrogen metabolism for metagenomic data using KofamScan(https://www.genome.jp/tools/kofamkoala/) ，NCyc and eggnog-mapper. But these results are significantly different. For example, the ORFs that are annotatated as gene nirK using three methods are different. Why is it so?

Thank you!

Dear @qichao1984

Thanks for the NCyc. This is me asking for guidance.

I'm using your script as follows:
$ perl NCycProfiler.pl -d ~/sickle/. -m diamond -f fastq -s nucl -si ~/sickle/samples.tsv -o ~/nitro_genes.tsv

After running, i get the following:

[...] Loading reference sequences... [1.2e-05s]
Deallocating buffers... [0.014192s]
Deallocating queries... [0.028869s]
Loading query sequences... [1e-05s]
Closing the input file... [7e-06s]
Closing the output file... [4.3e-05s]
Closing the database file... [5e-06s]
Deallocating taxonomy... [3e-06s]
Total time = 150.442s
Reported 153742 pairwise alignments, 153742 HSPs.
153742 queries aligned.
RP_IV_2 was not found in ~/sickle/samples.tsv, please check!

I tried running just two samples (two file names with their seqs counts in samples.tsv), though I got the same ¡error!

Information about samples.tsv

--> I have 10 fastq files in ~/sickle/. Those fastq files are listed in samples.tsv without the file extension, and their seq counts each. Values are \t separated.

So, Would you mind give a hand?

Thanks.
MaryoHg.

I am running the script as follows:

perl /programmes/NCyc/NCycProfiler.pl c -d /prepost/final_assemblies/matriline_assembly/contig_500bp_filt/cd-hit_dereplication_output/metaProdigal_output -m diamond -f faa -s prot -si seqnum.txt -o /prepost/final_assemblies/matriline_assembly/contig_500bp_filt/cd-hit_dereplication_output/metaProdigal_output/ncyc_output/ncyc_profiles.txt

yet I am getting no output to ncyc_profiles.txt. The run is executing correctly it seems and I do get the raw diamond output.

Any idea what I am doing wrong?

Hi,
I am interested in using NcycDB in my research but I see the Github repo does not contain a license for use of this code. I worry that means I cannot use the code in my research, because Github's support page says that no license means the code cannot be freely distributed.

Can you attach a license to this repo that allows for use of the code?

thanks for your time!

Hi, I attempted to run NCycProfiler.PL as follows:

$ perl NCycProfiler.PL -d ~/data/ -m /usr/bin/usearch -f fq -s nucl -si sample_info.txt -o NCyc_output
No diamond/usearch/blast files were detected!

What usearch files is the script trying to detect? I thought the script is intended to run usearch against the provided *faa.gz files. Am I supposed to run usearch (or diamond or blast) on my selected *faa.gz file before running the script, ie is the script looking for the output of a usearch run?

Here are some details on my usearch version, sample info file and read files

$ /usr/bin/usearch10
usearch v10.0.240_i86linux64, 65.6Gb RAM, 32 cores
(C) Copyright 2013-17 Robert C. Edgar, all rights reserved.                            http://drive5.com/usearch

License: [redacted], non-profit use, max 1 process(es)

$ head -2 sample_info.txt
17      2608559
18      325473

$ ls data/ | head -2
17.fq
18.fq

EDIT: updated usearch location and version info

I am getting output files looking like this but have no idea what is what?

k119_1934631_1  bur:Bcep18194_A4430     30.5    174     112     4       2       169     176     346     1.5e-12 68.2
k119_1254898_1  693661.Arcve_0699       32.6    215     121     4       33      223     14      228     1.1e-18 90.5
k119_1150324_1  5911.EAS00354   39.8    472     240     8       47      510     34      469     7.3e-89 282.7
k119_575171_1   1001530.PMSV_1997       46.8    111     51      1       1       103     1       111     2.4e-28 104.0
k119_1882347_1  34dd61e2fef5739d26ce3de0c76f9d94        48.5    171     86      2       1       169     176     346     1.8e-49 168.7
k119_1150325_2  6d0871a2edb457b604224c7994886b06        43.2    111     60      3       3       112     460     568     3.0e-16 77.8
k119_1516336_1  60051142d5b017a88152f3baddc63d81        36.6    243     127     6       29      244     238     480     1.4e-34 134.0
k119_1359477_1  374308263       71.4    56      16      0       1       56      1       56      1.1e-11 55.5
k119_1673198_1  350688.Clos_0484        68.9    405     120     4       4       405     383     784     9.9e-181        538.5
k119_1464047_1  339323450       49.5    107     53      1       235     340     42      148     7.9e-26 103.2
k119_522884_1   431895.XP_001747511.1   32.8    299     178     6       76      352     642     939     1.5e-42 160.6
k119_575174_1   866774.HMPREF9248_0150  78.1    32      7       0       9       40      23      54      3.0e-10 54.7
k119_1830061_1  642492.Clole_3012       66.2    287     97      0       1       287     387     673     2.4e-133        399.4
k119_1411762_1  212717.CTC00060 52.6    578     249     6       1       562     188     756     8.1e-202        585.5
k119_418314_1   e29be0bd0da4c63f6c0a39033ee3db32        48.9    47      24      0       5       51      93      139     3.4e-09 52.8
k119_470601_1   5825.PCHAS_101490       28.0    157     80      9       5       131     12      165     3.9e-06 49.3
k119_1359482_2  526221.XP_003003854.1   30.6    121     79      3       2       121     21      137     1.3e-08 53.5
k119_1202611_2  116747712       24.1    261     158     8       374     614     591     831     1.1e-05 52.0
k119_1202611_3  411464.DESPIG_00203     24.9    706     349     21      3       590     2       644     3.5e-42 165.2
k119_418316_2   a2fc249a6799579a69695a177760985a        53.0    66      31      0       168     233     10      75      3.5e-15 70.9
k119_366024_2   15669255        40.3    139     81      1       3       139     243     381     1.1e-26 105.1
k119_1307192_3  d9870872203d1ba76e85b283d4febb26        36.5    115     62      3       76      190     138     241     3.5e-08 54.7
k119_627464_1   fee9685fb9783d0114adb8c48bcfbf19        50.0    58      29      0       144     201     336     393     1.9e-06 50.1
k119_1254909_1  B6GAA7  62.3    154     58      0       1       154     461     614     1.1e-69 223.4
k119_418318_1   1ba5484590e9b2f67f7d6547226cc24c        48.1    131     67      1       1       130     97      227     3.4e-39 133.7

The first column is contig name and then I don't know what the rest are. Would you be able to provide some insight please?

Hi Professor Tu,

I set up blast using conda,

in the code,
##NCBI blast ftp://ftp.ncbi.nlm.nih.gov/blast/executables/legacy.NOTSUPPORTED/2.2.26/blast-2.2.26-x64-linux.tar.gz
(cannot get this gz file using this web)

my $blast = "~/bin/blast/bin/blastall";
(no blastall in bin folder, just have blastx, blastn, blastp in /home/lbhuang/miniconda3/pkgs/blast-2.12.0-pl5262h3289130_0/bin/)

my $formatdb = "~/bin/blast/bin/formatdb";
(no formatdb in bin folder)

so which one can I use?
Thank you very much

It is possible to get the classification for some interesting functional genes:

I want to know both the abundance and the classification for AOA, AOB, nirK/S and nosZ; wonder how can I do this using this database and pipeline.

thank you

Pretty much as title says. What is th enumber of sequences for the -si file?

Secondly, do you have any suggestions of best practices for the -rs options?

Hello,
I ran the perl script. The script runs till the following message is displayed:

Reported 91 pairwise alignments, 100 HSSPs
91 queries aligned.
#3
and terminates

I can't locate my output file. Can you mention where that might be.

Here is my script.

perl NCycProfiler.PL -d ~/prokka/cb2_bin9/ -m diamond -f fa -s nucl -si ~/prokka/cb2_bin9/ncyctry -o out_ncyc

Hi. :)
I am trying to run the following:

perl NCycProfiler.PL -d ~/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces -m diamond -f fasta -s prot -si /home/christian/Bioinformatic_programs/NCyc-master/FileForNcyc.txt -o PC_Ncyc

However, I get the following error:

cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10786.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10787.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10788.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10789.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10790.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10791.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10792.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10793.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10794.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10795.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10796.diamond': No such file or directory
cp: cannot stat '/home/christian/Desktop/scripts_to_abacus/From_abacus/PC_AA_Seqeunces/AAsequencesF10797.diamond': No such file or directory
was not found in /home/christian/Bioinformatic_programs/NCyc-master/FileForNcyc.txt, please check!

My tab delimited file looks like this:
head FileForNcyc.txt
AAsequencesF10786 18190
AAsequencesF10787 10828
AAsequencesF10788 41144
AAsequencesF10789 44519
AAsequencesF10790 874
AAsequencesF10791 81255
AAsequencesF10792 61666
AAsequencesF10793 79880
AAsequencesF10794 63433
AAsequencesF10795 51039

Maybe, the problem is trivial, however I cant see how to fix it? :)
Also, I tried with usearch, got the same problem

Cheers

I notice all the tools available here are alignment based. Would it not improve accuracy to use a HMMER/HMMSEARCH tool give this is metagenomic?

Secondly, might I suggest differentiating the working directory from the input directory? Essential add an option to direct where the NCyc database is stored?

I have used this on my metagenomic data but also want a profile for the host. Is this tool usable on eukaryotic genomes?

Hi there,

Could you please walk me through how to use this tool? I'm working on Mac Catalina and just installed perl.

Should I copy/paste the perl script into my terminal? Do I need to pre-download Diamond?

Many thanks!

Hello, teacher, I happened to see this tool developed by you, which is similar to the direction I want to do, but what I do is nitrogen metabolism of plants. I would like to ask, is this tool and data set suitable for my plant research?
Thanks for your reply.

Dear @qichao1984,

Thanks heaps for this DB.

We observed the following when analyzing our samples.
We have 10 assemblies (*.fa files) with different numbers of contigs.
I ran the NcyPerl.pl on my 10 assemblies (using the below command) and the obtained results contained a single hit with amoA_B in one sample (Figure 1).

$ NCycProfiler.pl -d ./ -m diamond -f fa -s nucl -si si.txt -rs 33175 -o out_file 2>&1 | tee log.run1

Figure 1: output file using all samples (n = 10) at the same depth (random sampling to 33,175)

However, after running NCycProfiler.pl on all my assemblies concatenated (sum of all contigs, around half a million), that hit is gone (Fig. 2). I cat my *.fa and run the script as shown below.

$ cat $(ls *.fa) >> cat_assembly.fa
$ NCycProfiler.pl -d ./ -m diamond -f fa -s nucl -si si.txt -o out_file 2>&1 | tee log.cat

Fig. 2. Output file using all samples (n = 10) concatenated (sum of all contigs)

We appreciate your time and response to this.
Greatings,

MaryoHg.

Cause I need to use the DNA seq of the N-related proteins, Can you publish the DNA sequence of each protein sequence?

Hello,

What are the units of the second column in the output file? Is this the number of sequences that mapped to that particular sequence? And, what is NR?
See image attached. Thanks!

I write your codes as following:
perl NCycProfiler.PL -d /home/lbertrans/ -m diamond -f fastq -s nucl -si Sequences_samples -o N_genes

The error is:
Please specify the location of diamond! at NCycProfiler.PL line 42.

I changed the path line 11:
my $diamond = "~/home/spack/opt/spack/linux-centos7-broadwell/gcc-4.9.4/diamond-2.0.11-obtc4vtc4dbd7irps2tz2jxstc4icaq5/bin/diamond"

It didn't work. I had the same error as before.

We are working in cluster server. Could it be the access to diamond software?

Thanks.

Lluís

Hello,
Thanks for providing this tool! I noticed that there are a number of proteins in NCyc_100_2019Jul.faa missing from id2gene.map.2019Jul (273,501 proteins vs. 219,091 entries in mapping file). Is this intentional?

Jeff

Hello!

According to the number of the second column in the output of the "NCycProfiler.PL" script, not all the annotation informations were found in the mapping file. How did this happen? (Software: usearch, database: ncyc_100_2019jul.faa, mapping file: id2gene.map.2019jul).

Thanks!

I read the NCycProfiler.PL, it seems that can align gene sequences to NCycDB database. But the sseqid (Subject Seqid) result of diamond blastx is not the protein ID in NCycDB (e.g. A0A254MU54). The sseqid (Subject Seqid) result of diamond blastx is a long string of values (e.g. 312309.VF). Is there any step I'm missing or how do I match the sseqid (Subject Seqid) of diamond blastx to the protein ID of NCycDB. Thanks!

Hi!
Thanks for creating this useful tool.
I have an error that I cant solve, even using the new script uploaded 20 days ago.

My command is
perl NCycProfiler2.PL -d . -m diamond -f .faa -s prot -si sampleinfi.txt -o output14

and apparenty eveything tus fine but this is the last lines of the file:

Total time = 917.505s
Reported 698857 pairwise alignments, 698891 HSPs.
698857 queries aligned.
cp: './proteins_aa14diamond' and './proteins_aa14diamond' are the same file
was not found in sampleinfi.txt, please check!

this is my -si
sampleinfi.txt

and the files I have in the folder with the script file are:
data NCycProfiler2.PL proteins_aa14.faa
diamond NCycProfiler.PL README.md
diamond_manual.pdf LICENSE proteins_aa14diamond sampleinfi.txt

Can you please, please help me?

Thanks in advance,

Ed

Hi! While the using of NCycDB, I found some potential inaccurate issues in it. So I want to confirm them here from the author for further improvement of this useful database.

1. There are 441 duplicated IDs (and sequence) in the file NCyc_100.faa. For your information, ten duplicated IDs are listed below:
   1000565.METUNv1_00600
   100226.SCO5525
   1002339.HMPREF9373_0600
   1002339.HMPREF9373_0712
   1003195.SCAT_4357
   1004836.SXCC_03526
   1006006.Mcup_0918
   1006551.KOX_23300
   1006551.KOX_27050
   1026882.MAMP_00572

2. There are 81 duplicated IDs in the file id2gene.map, and what's more, the classification of anyone of these IDs are not exactly the same. Ten of the mentioned IDs are listed below for your information.
   A0A060I505 napA
   A0A060I505 nasA
   A0A0D6FVV7 narG
   A0A0D6FVV7 narZ
   A0A127KYU8 nasB
   A0A127KYU8 nirB
   A0A142WVR2 narB
   A0A142WVR2 nasA
   A0A142Y6U0 narB
   A0A142Y6U0 nasA
   A0A142YH98 narB
   A0A142YH98 nasA
   A0A145CM87 narB
   A0A145CM87 nasA
   A0A1D8AZ75 narB
   A0A1D8AZ75 nasA

3. The ID list should be identical between the two files NCyc_100.faa and id2gene.map. But 54410 IDs in file NCyc_100.faa can not found in file id2gene.map, and 99288 IDs in file id2gene.map can not found in file NCyc_100.faa. The issue will cause inaccuracy when calculating the relative abundance of Nitrogen cycling genes (sub)families.

Should short reads from metagenomic data be mapped to those genes (or protein) (>90% identity for blastx or diamond blastx for example) that were found in the database to calculate relative abundance of each gene? but not only how many hits of that gene you have against the database because some shot reads may not be assembled successfully but they are actually part of those genes. I do not see any discussion about this in the SBB paper.