以Sample1-6的6个样品为例,二代双末端测序
(Sample1-fq_1.gz Sample1-fq_2.gz ... Sample6-fq_1.gz Sample6-fq_2.gz)
一 Row data质控
fastp -i Sample1-fq_1.gz -I Sample1-fq_2.gz -o Sample1-fq_1.clean.fq.gz -O Sample1-fq_2.clean.fq.gz -z 4 -q 20 -u 30 -n 10 -h Sample1.clean.html
(Sample2 - Sample6以同样的方式进行质控)
二 Hisat2比对获得bam文件
#单末端
hisat2 -p 2 -x Index -U read_1 | samtools sort -@ 20 > Sample1.sorted.bam
#双末端
hisat2 -p 2 -x index -1 read_1 -2 read_2 | samtools sort -@ 20 > Sample1.sorted.bam
三 利用FeatureCounts计算reads数
featureCounts -p -a gtf -g gene_id -o count.txt Sample1.bam Sample2.bam Sample3.bam Sample4.bam Sample5.bam Sample6.bam
# -g参数 根据gtf文件中基因ID的前缀
cut -f 1,7- count.txt | grep -v '#' > CountMatrix.csv
四 利用R进行定量分析(建议使用Rstudio-server)
library('DESeq2')
countdata <- read.table('CountMatrix.csv', row.names = 1,stringsAsFactors = T,check.names = F)
#CountMatrix.csv文件左上角为空
head(countdata)
coldata <- read.table('sample_table.txt',row.names = 1,stringsAsFactors = T)
#sample_table.txt 里面的顺序和CountMatrix.csv保持一致,注意condation里面是三个名一样,不是-1/-2/-3
head(coldata)
dds <- DESeqDataSetFromMatrix(countData = countdata,colData = coldata,design = ~ Condition)
dds <- DESeq(dds)
res <- results(dds)
head(res)
# 将res结果(差异倍数)与dds的counts(具体表达量)合并,并且按照行名进行合并,如果不加by="row.names",就不会按行名合并,那么就是如果6行文件和6行文件合并,最终结果为36行
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by="row.names",sort=FALSE)
nrow(resdata)
write.csv(resdata,file = "DEseq2_all_normalized.csv",row.names = FALSE)
###筛选差异基因###
table(res$padj<0.05) #统计padj<0.05的行数
table(res$log2FoldChange > 1 | res$log2FoldChange < -1) #统计log2FC >1 或 log2FC <-1的行数
res <- res[order(res$padj),]
diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #按条件进行筛选
resdata <- merge(as.data.frame(diff_gene_deseq2),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
nrow(resdata)
write.csv(resdata,file= "Rnaseq.diff.csv",row.names = F)
最终获得的Rnaseq.diff.csv包含了每个差异基因在各个样品中的表达量以及差异倍数