This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Visualize nucleotides for every position in R

Hi, please, how can I visualize nucelotides for every position. I have a percentage of bases but now I don't know how to create a plot.

fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
fastq
freq <- alphabetFrequency(fastq, as.prob = T,baseOnly=T)

I would like something like this (sorry for quality):

r
• 5,664 views • link
updated by xiaoleiusc xiaoleiusc • written by marija • marija

3 answers

You need to change your code. It is not alphabet frequency. It is consensus matrix. Example code and image below:

library(Biostrings)
fastq <- readDNAStringSet("sample.fastq","fastq")
## At each position, base frequency
afmc=consensusMatrix(fastq, baseOnly=T,as.prob = T)
tafmc=t(afmc)
matplot(tafmc[,-5], type="l", lwd=2, xlab="Read Length", ylab= "Base frequency at each position")
legend(legend = colnames(tafmc)[-5],"topright",col=1:4, lty=1:4, lwd=2)

test

5th column is others (other than A,T, G and C).

• 0 views • link
written by cpad0112 cpad0112

Awesome, thank you very much.

• 0 views • link
written by marija • marija

much simpler is using SeqTools in R:

library(seqTools)
# Reads fastq file
fq=fastqq("sample.fastq")
# Plots nucleotide frequency
plotNucFreq(fq,1)

Rplot

• 0 views • link
written by cpad0112 cpad0112

Thanks, does this apply to fasta file containing sequence of one gene? How to set sliding window size in Seqtools?

• 0 views • link
written by xiaoleiusc xiaoleiusc 
library(Biostrings)
fastq <- readDNAStringSet("https://d28rh4a8wq0iu5.cloudfront.net/ads1/data/ERR037900_1.first1000.fastq","fastq")
freq <- alphabetFrequency(fastq, as.prob = T,baseOnly=T)
matplot(freq,type='l')
legend(legend = colnames(freq),x=700,y=0.7,lty=1:5,col=1:5)

Note answer about consensusMatrix. That code is the best answer.

• 0 views • link
written by Sean Davis Sean Davis

This is per-cycle base content. R is not good at, and is not designed for such task.

Why complicate life? Just use fastp to do the easiest and fastest way for such tasks, including QC and filtering. It will output the figures like what you want.

• 0 views • link
written by chen chen

Log in to answer this question.

More posts like this