miRNA-seq分析流程

轉(zhuǎn)載2016-11-28 16:29:44

miRNA是生物中非常重要的一類非編碼小RNA，其在生物體的調(diào)控中具有非常重要的作用，在人中大約三分之一的基因受到miRNA的調(diào)控。對(duì)于miRNA轉(zhuǎn)錄后調(diào)控的分析也越來(lái)越多。那么拿到一組miRNA測(cè)序的數(shù)據(jù)之后我們進(jìn)行怎樣的分析呢？

第一，??對(duì)于所有的測(cè)序數(shù)據(jù)，我們都要進(jìn)行質(zhì)量的檢測(cè)，這里我常用的檢測(cè)軟件是fastQC，（fastqcdata1.fq -o data1），得到的結(jié)果是一個(gè)文件夾的壓縮形式，里面可以得到以下所示的信息：

網(wǎng)頁(yè)結(jié)果展示

?同時(shí)這些結(jié)果都有單獨(dú)的圖片格式，用于數(shù)據(jù)展示與質(zhì)量評(píng)定。在新版本的fastQC中有一個(gè)新的功能就是識(shí)別reads中包含的adapter序列，并且fastqQC中有一個(gè)adapter的庫(kù)，可以從里面找到對(duì)應(yīng)的adapter序列，再也不用擔(dān)心沒有測(cè)序報(bào)告沒法去掉adapter了。

第二，??對(duì)數(shù)據(jù)進(jìn)行過(guò)濾，這里我用的是cutadapt，這個(gè)軟件可以去掉reads中的adapter，低質(zhì)量的reads以及過(guò)長(zhǎng)過(guò)短的reads，還可以對(duì)reads中含有N的進(jì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，這里對(duì)clean的reads還要進(jìn)行一下長(zhǎng)度分布的統(tǒng)計(jì)，一般就是自己寫腳本，我用的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)

統(tǒng)計(jì)完長(zhǎng)度分布之后就是做呈現(xiàn)出來(lái)，這里是用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的位置有兩個(gè)峰：

第四，??在得到高質(zhì)量的clean數(shù)據(jù)之后就是進(jìn)行比對(duì)，將miRNA的數(shù)據(jù)比對(duì)到相應(yīng)物種的基因組上，這里我用的是bowtie軟件，（bowtie -q -v 2 -l 10 -k 15 Reference/genome.fa trim_data1.fq -S data1.sam 2>mapping.info），我分析的植物miRNA-seq的數(shù)據(jù)，比對(duì)率超過(guò)了90%。

第五，??在得到比對(duì)的結(jié)果之后就是用HTSeq進(jìn)行count計(jì)數(shù)，把在不同的材料中表達(dá)的miRNA的reads支持?jǐn)?shù)統(tǒng)計(jì)出來(lái)。

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

第六，??然后就是重頭戲，差異表達(dá)的miRNA，這里分為有重復(fù)的處理和沒有重復(fù)的處理兩種，對(duì)于沒有重復(fù)的用DEGseq處理，有重復(fù)的用DESeq處理。

沒有重復(fù)的用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()

?有重復(fù)的用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()

第七，??在找到相應(yīng)的差異表達(dá)miRNA之后，對(duì)其靶基因進(jìn)行預(yù)測(cè)與分析，這里我是將我做的玉米miRNA的所有靶基因進(jìn)行了預(yù)測(cè)，這里的玉米miRNA的全部注釋信息是從miRbase中下載的，得到miRNA靶基因的詳細(xì)信息之后，對(duì)于不同組的差異表達(dá)miRNA可以從中去對(duì)應(yīng)分析。

在植物中有兩個(gè)比較好的miRNA預(yù)測(cè)的軟件，分別是psRNAtarget和psRobot，psRNAtarget是一個(gè)支持在線預(yù)測(cè)的軟件，psRobot可以將軟件安裝在服務(wù)器中，用命令行進(jìn)行預(yù)測(cè)。

在得到的結(jié)果中兩種軟件預(yù)測(cè)到的共同的部分是結(jié)果比較可信的部分。

第八，??當(dāng)然對(duì)于miRNA還有很過(guò)方面，對(duì)靶基因的功能分析，對(duì)miRNA二級(jí)結(jié)構(gòu)的分析，對(duì)樣品中新miRNA的分析等等。