This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Trimming an alignment in R

Hi! I was wondering if anyone could help. I'm using the package ape in R and I want to trim an alignment of class DNAbin, this is, selecting portions of the alignment between two nucleotide positions. For example, I want all the sequences of my alignment but only between the nt positions 1503 and 3051.

Is it possible with ape or any other R package?

Thanks!

G.

ape r alignment
• 7,204 views • link
updated by Ram Ram • written by Gon • Gon

Thanks Andrew! Eventually I managed to do it with ape, converting the object type DNAbin in a matrix, that way you can subset the usual way with square brackets. If you initially import the alignment from a fasta file using the argument as.matrix=TRUE you can skip the step of converting the DNAbin to a matrix. In my example:

full_alignment <- read.dna("MyFastaFile.fas", format="fasta", as.matrix=TRUE)
gene_wanted <- full_alignment[,1503:3051]
• 0 views • link
updated by Ram Ram • written by Gon • Gon

Hi Gon, how did you then write the subset to a new fasta file, without losing the fasta structure?

• 0 views • link
written by jbt38 • jbt38

Hi jbt38! Following my example, I would export the sliced alignment into a new file like this:

write.dna(gene_wanted, file="gene_wanted.fas", format="fasta", nbcol=-1, colsep="")

Hope it's useful.

• 0 views • link
written by Gon • Gon

Legend, thanks for this - especially remembering what you did from so long ago!

• 0 views • link
written by jbt38 • jbt38

1 answer

I'm pretty sure that all you want to do is subset your alignment right? You can use samtools view:

samtools view -h /path/to/file.bam 17:7512445-7513455
• 0 views • link
written by andrew.j.skelton73 andrew.j.skelton73

Log in to answer this question.

More posts like this