Thanks very much. I will get all unique names first then.
• 0 views
•
link
written
by
phxu6780 •
duplicate gene IDs with different length after cd-hit-est dereplicate and cluster
Hi. I am trying to use CD-hit to remove the duplicates from the file that is the combine of all nucleic acid fastq output from prodigal. I used the following parameters:
cd-hit-est -i nuc_sum.fa -o cd-hit_sum -c 0.95 -s 0.8 -M 0 -T 0 -n 8
The representative sequencing shows in the fasta file. But there are many small or fragmented sequence with the same header ID.
Does anyone know how to set the parameters in cd-hit-est to make sure there will be only one sequence for one header ID?
I tried cd-hit-est -i nuc_sum.fa -o cd-hit_sum -c 1 -t 1 -d 0, which someone recommend. But it does not solve this problem.
• 6,795 views
•
link
written
by
phxu6780 •
2 answers
CD-HIT removes sequence redundancy, not ID redundancy. So if there are repeated IDs after clustering, that is because you have multiple IDs that actually contain different sequences, at least given the clustering parameters.
First, I suggest you used the -d 0 switch, which will not truncate fasta header descriptions. That way you can be sure that the names are truly identical. If multiple names still persist, i don't think you will get CD-HIT to fix that because it doesn't look into IDs at all. If you start with all unique names before clustering, you will not have this problem.
• 0 views
•
link
written
by
Mensur Dlakic
Do you have any suggestions to get unique names before clustering?
• 0 views
•
link
written
by
phxu6780 •
Great. If possible, I would like to keep all duplicate sequencing. I will try this suggestion first and see how the output will be in the downstream analysis.
• 0 views
•
link
written
by
phxu6780 •
There are two other ways to handle it:
Use
seqkit rmdup -n. This would remove sequences with identical names. However, I don't know how it works once it finds identical IDs (first encounter or retain one with longest sequence). Then feed the resultant fasta to CD-HIT.Use
awkorseqkitto convert fasta to tab file, group by column 1 (ID) and then pick the longest one in 2nd column, from each group, then convert back to fasta. Then feed the fasta to CD-HIT.
Sounds promising way to get unique IDs before cd-hit. Do you know the exact command of awk or seqkit to run this?
• 0 views
•
link
written
by
phxu6780 •
for seqkit:
$ seqkit rmdup -ni <input.fa> -o <output.fa>
If you want to save duplicated sequences in a separate file:
$ seqkit rmdup -ni -o <new.fa> -d <duplicated.fa> -D <duplicated_details.txt>
You need to use new.fa for further analysis. Duplicated.fa will give you duplicated sequences. This is for your reference only. seqkit retains first encounter, not the longest one.
I will post the awk code after some time.
Awesome. Thanks very much. I would like to get the longest one. Looking forwards to your post then.
• 0 views
•
link
written
by
phxu6780 •
with datamash and awk (datamash is available in most of the debian,*buntu, brew, conda repos):
input data
$ cat test.fa
>seq1
atgc
catg
>seq1
gtac
cagtgacg
>seq1
atgcgctgacgacgagcg
>seq11
atgc
>seq11
gatc
output:
$ awk '{if(NR==1) {print $0} else {if($0 ~ /^>/) {print "\n"$0} else {printf $0}}} END {print "\r"}' test.fa | awk -v OFS="\t" -v RS=">" 'NR>1 {print $1,$2, length($2)}' | datamash -sg1 max 3 -f | awk -v OFS="\n" '{print ">"$1,$2}'
>seq1
atgcgctgacgacgagcg
>seq11
gatc
with awk only (not tested):
$ awk '{if(NR==1) {print $0} else {if($0 ~ /^>/) {print "\n"$0} else {printf $0}}} END {print "\r"}' test.fa| awk -v OFS="\n" -v RS=">" ' NR>1 {if (length(a[$1]) < length($2))a[$1]=$2;}END{for(i in a){print ">"i,a[i];}}'
If sequences in fasta files are in single line, you can remove '{if(NR==1) {print $0} else {if($0 ~ /^>/) {print "\n"$0} else {printf $0}}} END {print "\r"}' part of the code. This would linearize the fasta file.
Log in to answer this question.
More posts like this
-
Why cd-hit-est not work when sequence identity threshold<0.95？
written by JZX •
-
Word size CD-HIT-EST
written by Nathan •
-
understanding CD hit command
written by avinashbelure3 •
-
PBS scripts for CD-HIT
written by lucia •
-
how to speed up cdhit clustering?
written by bitpir
-
How can run cd-hit-est with a clstr threshold less than 0.8?
written by m.koohi.m
-
how to use CD_HIT to remove the redundant sequence from trinity output file
written by wu.zhiqiang.1020 •
-
Trouble installing cd-hit
written by nlehmann
-
Trinity assembly and ORF decoder issue
written by onspotproductions
-
Cd-Hit-Est-2D Segmentation Fault
written by Mariam Rizkallah •
Please don't post data/code examples as images. Post actual data and format it using
101button in editor. It makes it difficult for people to grab the example (not in this case) to provide suggestions.Thanks very much for your suggestion, it is my first time to post questions. I will try to use editor another time.
Are you certain that is the case? Perhaps there is more stuff beyond
_1than what we see in image?I am sure it is the case. As I used bowtie2 to align clean raw reads to the index made by this fasta file cd-hit_sum. And I got cd-hit_sum.sam file. But when I ran
samtools view -b -S samfile.sam -o bamfile.bam, it turned out fail as reason of the duplicate entry in sam header. Then I figure out the reason should be that the cd-hit_sum.fa has many sequences with the same header, but the length of them are different.Check your
CD-HIToutputFASTAfile, and check whether you have more than one representative sequence from each of your headers. So something likegrep -oP "^>k[0-9]+_[0-9]+_[0-9]+" cd-hit_sum | sort | uniq -c.The information depicted in your screenshot here is the cluster output, and it just indicates what sequences belong to each cluster.
Hi. Thanks for suggestions. I tried
grep -oP 'k141_158430_1' cd-hit_sum | sort | uniq -c, the output is 9 k141_158430_1.So you have nine sequences that share the identifier
k141_158430_1(or at least that particular portion of it) that aren't/weren't similar enough to be clustered away under the settings you calledCD-HITwith. As @Mensur Dlakic mentioned,CD-HITdoes not cluster on the basis of sequence identifiers but on the basis of the sequences themselves.Now I am running again of cd-hit-est adding the command of -G 0 -aS 0.9.
-aS, alignment coverage for the shorter sequence, default 0.0 if set to 0.9, the alignment must covers 90% of the sequence.
Hope it will help. But not sure is this the case to make sure the output will have no same header for each gene or cluster.
Could you perhaps elaborate on the provenance of your data? Why is it that you wish to have one representative sequence per identifier? If one representative sequence per identifier that you strictly want--and want some clustering on top of that--I suggest you group your sequences by whatever happens to be the common identifier string, and select the longest sequence as the representative (or something to that effect), and follow that up with a light clustering step.