Trimming an alignment in R

Question

[Gon]

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.

[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]

[jbt38]

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

[Gon]

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.

[jbt38]

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

Answer 1

[andrew.j.skelton73]

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