这篇是Hisat2+HTSeq+DESeq2的流程。
首先补充一个说明，stringtie提供了一个叫prepDE.py的脚本，可以用stringtie的结果输出DESeq2需要的矩阵。 在rna-seq的第一篇中已经说过怎么下载了。 使用的方法是，先创建一个列表，列表形式点击这里查看。
# -g表示输出基因结果，-t表示输出转录本结果
python prepDE.py \
-i sample_list.txt \
-g gene_results.csv \
-t transcript_results.csv
下面是利用htseq-count来统计Hisat2比对之后的结果：
# 有多少个sample就写多少个
# -i 表示输出的是chrX.gtf中的哪个值
htseq-count -f bam -r pos -s no \
-i gene_id \
sample1.sort.bam \
sample2.sort.bam \
sample3.sort.bam \
sample4.sort.bam \
chrX.gtf \
1>sample_chrX.count \
2>sample_htseq.log
然后我用脚本给输出的文件加上头和删掉后面的信息行：
#!/usr/bin/python2.7
countfile = open("sample_chrX.count", "r")
results = open("sample_chrX_results.csv", "w")
lines = countfile.readlines()
aft = lines[:-5]
head = "这里填入要加进去的标题，用tab分割"
aft.insert(0, head)
reuslts.writelines(aft)
reuslts.close()
countfile.close()
最后就来到DESeq2的操作步骤，使用R：
library(DESeq2)
# 数据预处理
database <- read.table(file = "sample_chrX_results.csv", sep = ",", header = TRUE, row.names = 1)
database <- round(as.matrix(database))
# 设置分组信息并构建dds对象，按照geuvadis_phenodata.csv给的信息写
condition <- factor(
c(
"male",
"male",
"female",
"female",
"male",
"female",
"female",
"male",
"male",
"female",
"male",
"female"
),
levels = c("male", "female")
)
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
dds <- dds[ rowSums(counts(dds)) > 1, ]
# 使用DESeq函数估计离散度，然后差异分析获得res对象
dds <- DESeq(dds)
res <- results(dds)
# 最后设定阈值，筛选差异基因，导出数据(全部数据。包括标准化后的count数)
res <- res[order(res$padj),]
diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene <- row.names(diff_gene)
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)
write.csv(resdata,file = "male_female_chrX_diff.csv",row.names = FALSE)
作图还没学会！等学会了再来说怎么作图。