原來不用繁瑣得安裝調(diào)試Linux系統(tǒng)，在Windows下就能快速分析二代測序轉(zhuǎn)錄組數(shù)據(jù)。Rsubread就能實現(xiàn)此功能。說明書里面的圖很直觀，直接上圖了（圖來源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/