miRNA-seq分析流程
转载2016-11-28 16:29:44
miRNA是生物中非常重要的一类非编码小RNA，其在生物体的调控中具有非常重要的作用，在人中大约三分之一的基因受到miRNA的调控。对于miRNA转录后调控的分析也越来越多。那么拿到一组miRNA测序的数据之后我们进行怎样的分析呢？
第一， 对于所有的测序数据，我们都要进行质量的检测，这里我常用的检测软件是fastQC，（fastqcdata1.fq -o data1），得到的结果是一个文件夹的压缩形式，里面可以得到以下所示的信息：
网页结果展示
同时这些结果都有单独的图片格式，用于数据展示与质量评定。在新版本的fastQC中有一个新的功能就是识别reads中包含的adapter序列，并且fastqQC中有一个adapter的库，可以从里面找到对应的adapter序列，再也不用担心没有测序报告没法去掉adapter了。
第二， 对数据进行过滤，这里我用的是cutadapt，这个软件可以去掉reads中的adapter，低质量的reads以及过长过短的reads，还可以对reads中含有N的进行处理。（cutadapt-a AGATCGGAAGAG --quality-base 33 -m 10 -q 20 --discard-untrimmed -o trim_data1.fqdata1.fq > cutadpt.info），这里--discard-untrimmed是把reads中不含有adapter的reads去掉。
第三， 由于分析的是miRNA-seq，这里对clean的reads还要进行一下长度分布的统计，一般就是自己写脚本，我用的python。
importsys
miRNA_len= {}
fori in range(0,52):
miRNA_len[i] = 0
fori in open(sys.argv[1]):
if i.startswith('@') ori.startswith('+'):
continue
length = len(i) - 1
miRNA_len[length] += 1
fori in miRNA_len:
printstr(i)+"\t"+str(miRNA_len[i]/2)
统计完长度分布之后就是做呈现出来，这里是用R作图：
#use all reads
s1= read.table("trim_data1.stat")
s2= read.table("trim_data2.stat")$V2
s3= read.table("trim_data3.stat")$V2
s4= read.table("trim_data4.stat")$V2
s5= read.table("trim_data5.stat")$V2
s6= read.table("trim_data6.stat")$V2
data= cbind(s1, s2, s3, s4, s5, s6)
colnames(data)= c("Length","data1","data2","data3","data4","data5","data6")
#normalize by library size
data$kBT_0= 100 * data$data1/sum(data$data1)
data$kBT_1= 100 * data$data2/sum(data$data2)
data$kBT_3= 100 * data$data3/sum(data$data3)
data$kN6_0= 100 * data$data4/sum(data$data4)
data$kN6_1= 100 * data$data5/sum(data$data5)
data$kN6_3= 100 * data$data6/sum(data$data6)
library(reshape2)
data.melt= melt(data, id="Length")
library(ggplot2)
p<- ggplot(data.melt, aes(x=Length, y=value, col=variable))
p+ geom_line() +
theme( text = element_text(size=30),
panel.background=element_blank(),
axis.line = element_line(size = 1,colour="black"),
axis.text =element_text(colour="black")) +
labs(title="All reads lengthdistribution",x="Read Length", y="Fraction (%)")
得到的示意图如下，一般在21和24nt的位置有两个峰：
第四， 在得到高质量的clean数据之后就是进行比对，将miRNA的数据比对到相应物种的基因组上，这里我用的是bowtie软件，（bowtie -q -v 2 -l 10 -k 15 Reference/genome.fa trim_data1.fq -S data1.sam 2>mapping.info），我分析的植物miRNA-seq的数据，比对率超过了90%。
第五， 在得到比对的结果之后就是用HTSeq进行count计数，把在不同的材料中表达的miRNA的reads支持数统计出来。
for i in data1 data2 data3 data4 data5 data6
do
htseq-count-s no -t miRNA -i ID -o $i.hc.sam $i.ht.sam miRNA_reference/miRNA.gff3 | tee$i.count &
ls$i.count >> count.list
done
第六， 然后就是重头戏，差异表达的miRNA，这里分为有重复的处理和没有重复的处理两种，对于没有重复的用DEGseq处理，有重复的用DESeq处理。
没有重复的用DEGseq：
R
library("DEGseq")
#BT_0_1
geneExpMatrix1<- readGeneExp("ht.genotype_data1.txt", geneCol=1, valCol=3)
geneExpMatrix2<- readGeneExp("ht.genotype_data2.txt", geneCol=1, valCol=2)
write.table(geneExpMatrix1[30:31,],row.names=FALSE)
write.table(geneExpMatrix2[30:31,],row.names=FALSE)
pdf(file="data1_2.pdf")
layout(matrix(c(1,2,3,4,5,6),3, 2, byrow=TRUE))
par(mar=c(2,2, 2, 2))
DEGexp(geneExpMatrix1=geneExpMatrix1,geneCol1=1, expCol1=2, groupLabel1="data1",
geneExpMatrix2=geneExpMatrix2,geneCol2=1, expCol2=2, groupLabel2="data2",
method="MARS",outputDir="05DEmiRNA/DEGSeq")
dev.off()
有重复的用DESeq：
R
library("DESeq")
data=read.table("ht.genotype_data.txt",header=TRUE,row.names=1)
pd=data.frame(row.names=colnames(data),condition=c("data3","data3","data4","data4"),libType=c("single-end","single-end","single-end","single-end"))
ps=pd$libType=="single-end"
ct=data[,ps]
condition=pd$condition[ps]
cds=newCountDataSet(ct,condition)
cds= estimateSizeFactors(cds)
sizeFactors(cds)
cds= estimateDispersions(cds)
res=nbinomTest(cds,"data3","data4")
write.table(res,file="data3_data4.xls")
quit()
第七， 在找到相应的差异表达miRNA之后，对其靶基因进行预测与分析，这里我是将我做的玉米miRNA的所有靶基因进行了预测，这里的玉米miRNA的全部注释信息是从miRbase中下载的，得到miRNA靶基因的详细信息之后，对于不同组的差异表达miRNA可以从中去对应分析。
在植物中有两个比较好的miRNA预测的软件，分别是psRNAtarget和psRobot，psRNAtarget是一个支持在线预测的软件，psRobot可以将软件安装在服务器中，用命令行进行预测。
在得到的结果中两种软件预测到的共同的部分是结果比较可信的部分。
第八， 当然对于miRNA还有很过方面，对靶基因的功能分析，对miRNA二级结构的分析，对样品中新miRNA的分析等等。
