Windows下分析二代测序RNASeq数据快速入门

原来不用繁琐得安装调试Linux系统，在Windows下就能快速分析二代测序转录组数据。Rsubread就能实现此功能。说明书里面的图很直观，直接上图了（图来源https://bioconductor.org/packages/release/bioc/vignettes/Rsubread/inst/doc/SubreadUsersGuide.pdf）

代码如下：

setwd("G:/Rsubread_demo")

library(Rsubread)

buildindex(basename="ATG_index",reference="Arabidopsis_thaliana.TAIR10.dna_sm.chromosome.1.fa")

read1.files <- list.files(path = "E:/迅雷下载", pattern = ".+\\_1.fastq.gz", ignore.case = TRUE,

                        full.names = TRUE, recursive = TRUE) 

read2.files <- list.files(path = "E:/迅雷下载", pattern = ".+\\_2.fastq.gz", ignore.case = TRUE,

                          full.names = TRUE, recursive = TRUE) 

bamfiles=paste(stringi::stri_sub(read1.files,1,18),"_1.fastq.gz.subread.BAM",sep="")

align(index="ATG_index",readfile1=read1.files,readfile2=read2.files,nthreads=)

fc <- featureCounts(files=bamfiles,annot.ext="Arabidopsis_thaliana.TAIR10.36.gtf",isGTFAnnotationFile=T)

fc$stat

samplename <- colnames(fc$counts)

condition<- c("d0","d3","d3","d3","d3","d0","d0","d0")

samples=data.frame(samplename,condition)

samples

library(edgeR);library(limma)

x <- DGEList(fc$counts,lib.size = colSums(fc$counts),genes=fc$annotation[,c("GeneID","Length")],

             norm.factors = calcNormFactors(fc$counts),samples = samples$samplename,remove.zeros=T, group = samples$condition)

isexpr <- rowSums(cpm(x) > 10) >= 2

x <- cpm(x[isexpr,],log=T)

treat_time <- factor(samples$condition)

design <- model.matrix(~0+treat_time)

colnames(design) <- levels(treat_time)

rownames(design)=samples$samplename 

group <- factor(samples$condition)

# View an MDS plot

plotMDS(x, labels=treat_time, col=as.numeric(treat_time),xlab="Dim1",ylab="Dim2")

fit <- lmFit(x,design)

contr <- makeContrasts(d3-d0,levels=design)

fit.contr <- eBayes(contrasts.fit(fit,contr))

dt <- decideTests(fit.contr,p.value = 0.05,adjust.method = "none",lfc = 1.3)

summary(dt)

options(digits=2)

topTable(fit.contr)

volcanoplot(fit.contr, coef=1,style = "p-value", highlight = 100, names = fit.contr$GeneID,   hl.col="blue", xlab = "Log2 Fold Change")

注：

gtf下载地址：ftp://ftp.ensemblgenomes.org/pub/release-36/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.36.gtf.gz

fa下载地址：ftp://ftp.ensemblgenomes.org/pub/release-36/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna_sm.toplevel.fa.gz

RNASeq下载地址：https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-5130/samples/