1
$\begingroup$

I have a Multiple Sequence Alignment, where there are around 5000 sequences. There also exists many sequences where, there are so much of non-sequenced regions (for instance, AU----CGGGCA--NNNNNNNNNN).

So, by trial and error, I found a threshold value of 25%, so I removed the sequences which contained more than 25% of non-sequenced region (N's here), retaining around 400 sequences in the alignment, by writing a simple R script.

But, I would like to know, if there is any other standard method to do this, because I do not think, removing sequences based on finding a threshold for the percentage of non-sequenced regions, by trial and error is more efficient.If there is, can you please tell me, how to do?

Thanks in advance.

Improve this question
$\endgroup$
10
  • 2
    $\begingroup$ To remove sequences that contain ambiguous sequence is pretty straightforward. You can do it programmatically instead of doing manually. What do you mean by proper method? Can you please clarify what your objective is? $\endgroup$
    – WYSIWYG
    Commented Jul 25, 2019 at 8:43
  • $\begingroup$ Sure, I removed the sequences with a simple R script. I do not think, removing sequences based on finding a threshold for the percentage of non-sequenced regions, by trial and error is more efficient. I was wondering if there is any standard way to do it. $\endgroup$
    – User
    Commented Jul 25, 2019 at 8:49
  • $\begingroup$ That totally depends on your judgment. There can be sequence blocks that do not contain ambiguous bases. You can take them if you think the entire sequence is not necessary. This could be done only for local alignments. For doing global alignments, you should take a call on whether or not the number of ambiguous bases is too much. It would also matter if the ambiguous bases are present together or spread out. A stretch of Ns would be worse than the Ns being spread out. You can try MSA with and without these sequences and see if there is any significant change in the consensus. $\endgroup$
    – WYSIWYG
    Commented Jul 25, 2019 at 8:58
  • 3
    $\begingroup$ Deciding the threshold depends on your prior knowledge about the sequence region. The software is unaware of the biology. For e.g. the N could be in non-conserved regions: it would not probably affect the MSA. If they are in the conserved regions, they could distort the analysis. Anyway, can you please edit your question to clarify the exact objective (as you mentioned in comments)? Comments are not searchable and can go unnoticed. Moreover, the post should be self-sufficient. $\endgroup$
    – WYSIWYG
    Commented Jul 25, 2019 at 9:02
  • $\begingroup$ I edited my question. $\endgroup$
    – User
    Commented Jul 25, 2019 at 9:19

2 Answers 2

Reset to default
1
$\begingroup$

The bioconductor package "ShortRead" might be what you are after.

It has a function "clean" to remove sequences with ambiguous bases. It also can remove duplicates if you need to.

Improve this answer
$\endgroup$
1
$\begingroup$

It would be helpful to know your application. For some applications, having ambiguous sequences is not necessarily a problem as long as you have a mechanism to handle such non-informative positions. The information from the other positions may still be valuable.

An alternative to removing the sequences (rows of the alignment matrix) with many Ns would be to remove non-informative or poorly aligned positions (columns of the alignment matrix) from the alignment. It depends to some extent on what your goal is. Tools such as gblocks can be used in this way.

For example, if your goal is to perform clustering or phylogenetic estimation of the region, or identify highly conserved regions, then removing poorly aligned or non-informative columns would be worthwhile.

On the other hand, if your goal is somehow related to the entire sequence (e.g. protein structures or open reading frames), then removing rows/sequences may make more sense. But it is still not clear to me that you really need to remove ambiguous positions. For some discussion of this you can see this old researchgate thread or this biostars thread.

Improve this answer
$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .