How to put characters

Question

I have a fasta format file where in i have to only keep those nodes whose length is less than 100. however, the problem i am currently facing is that i am able to separate the nodes but am not able to put the characters of each node in separate variable whose length i can then check and subsequently separate the requisite nodes from longer ones. So what i mean is i am able to read the headings and separate nodes but how do i put the characters within each node in a variable.

This is a sample of my data

>NODE_1
GTTGGCCGAGCCCCAGGACGCGTGGTTGTTGAACCAGATCAGGTCCGGGCTCCACTGCAC
GTAGTCCTCGTTGGACAGCAGCGGGGCGTACGAGGCCAGCTTGACCACGTCGGCGTTGCG
CTCGAGCCGGTCATGAACGCGGCCTCGGCGAGGGCGTTCTTCCAGGCGTTGCCCTGGGAA

>NODE_2
CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGCGCACCGACCTGTCCA
CCGCGGACAAGGCCGGTTACCTGCACCGCTACATCGAGGCCAGCCGCATCGCGTTCGCGG
ACCGCGGGCGCTGGGTCGGCGACCCCGCCTTCGAGGACGTAC

>NODE_3
CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGCGCACCGACCTGTCCA
CCGCGGACAAGGCCGGTTACCTGCACCGCTACATCGAGGCCAGCCGCATCGCGTTCGCGG
ACCGCGGGCGCTGGGTCGGCGACCCCGCCTTCGAGGACGTACATCATTCCTTAATCTTCC

my code:

x <- readLines("1.fa", n = -1L, ok = TRUE, warn = TRUE)

for (i in 1:length(x)) {
    if (substr(x[i],1,1)=='>') {
        head <- c(head,x[i])
        q <- x[i+1] 
        if (q=!0) {
            contig <- c(contig,q)
            print(contig)       
            contig.length <- c(contig.length, nchar(q))
        } else {
            break
        }
    } else {
        z <- paste(z,x[i], sep=" ")
    }
}

Answer 1

You should use BioConductor for that. You're actually trying to parse a FASTA-file to some kind of a list. Bioconductor has a simple function read.fasta() that does just that, and returns an object where you can get the lengths and so on. Learning bioconductor is definitely worth the hassle if you work with sequences.

To do it in base R, you'll need to work with lists, something like :

Split.Fasta <- function(x){
  out <- list()
  for(i in x){
    if(substr(i,1,1)==">") {

      name <- gsub(">","",i)
      out[[name]] <- character(0)

    } else if (grepl("\\w",i)){
      out[[name]] <- paste(out[[name]],gsub("\\W","",i),sep="")
    } 
  } 
  out
}

Which works like :

zz <- textConnection(">NODE_1 
GTTGGCCGAGCCCCAGGACGCGTGGTTGTTGAACCAGATCAGGTCCGGGCTCCACTGCAC 
GTAGTCCTCGTTGGACAGCAGCGGGGCGTACGAGGCCAGCTTGACCACGTCGGCGTTGCG 
CTCGAGCCGGTCATGAACGCGGCCTCGGCGAGGGCGTTCTTCCAGGCGTTGCCCTGGGAA

>NODE_2 
CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGCGCACCGACCTGTCCA 
CCGCGGACAAGGCCGGTTACCTGCACCGCTACATCGAGGCCAGCCGCATCGCGTTCGCGG 
ACCGCGGGCGCTGGGTCGGCGACCCCGCCTTCGAGGACGTAC

>NODE_3 
CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGCGCACCGACCTGTCCA
CCGCGGACAAGGCCGGTTACCTGCACCGCTACATCGAGGCCAGCCGCATCGCGTTCGCGG
ACCGCGGGCGCTGGGTCGGCGACCCCGCCTTCGAGGACGTACATCATTCCTTAATCTTCC")

X <- readLines(zz,n=-1L,ok=TRUE,warn=TRUE)
close(zz)

Y <- Split.Fasta(X)
$`NODE_1 `
[1] "GTTGGCCGAGCCCCAGGACGCGTGGTTGTTGAACCAGATCA...

$`NODE_2 `
[1] "CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGC...

$`NODE_3 `
[1] "CCTCCGGCGGCACCACGGTCGGCGAGGCCCTCAACATCCTGGAGCGCAC...

It returns a list which you can use later on to check lengths and so on :

sapply(Y,nchar)
NODE_1  NODE_2  NODE_3  
    180     162     180

Still, learn to use BioConductor, you'll thank yourself for that.

Answer 2

You could install the seqinr package, which has lots of methods for analysing sequence data.

install.packages("seqinr")

Next, read in your fasta file:

seqs <- read.fasta("myfile.fa")

And then, extract sequences from the list with length < 100:

seqs.small <- seqs[sapply(seqs, function(x) getLength(x) < 100)]