miRNA-seq分析流程

最新推荐文章于 2024-05-23 09:39:52 发布

高锦-生信

最新推荐文章于 2024-05-23 09:39:52 发布

阅读量1.7w

点赞数 3

文章标签： miRNA-seq分析流程

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 inrange(0,52):

miRNA_len [i] = 0

fori inopen(sys.argv [1] ):

if i.startswith('@') ori.startswith('+'):

continue

length = len(i) - 1

miRNA_len[length] += 1

fori inmiRNA_len:

printstr(i)+"\t"+str(miRNA_len [i] /2)

统计完长度分布之后就是做呈现出来，这里是用R作图：

#use allreads

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")

#normalizeby 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 readslengthdistribution",x="Read Length", y="Fraction (%)")

得到的示意图如下，一般在21和24nt的位置有两个峰：

第四，在得到高质量的clean数据之后就是进行比对，将miRNA的数据比对到相应物种的基因组上，这里我用的是bowtie软件，（bowtie -q -v 2 -l 10 -k 15 Reference/genome.fa trim_data1.fq -Sdata1.sam 2>mapping.info），我分析的植物miRNA-seq的数据，比对率超过了90%。

第五，在得到比对的结果之后就是用HTSeq进行count计数，把在不同的材料中表达的miRNA的reads支持数统计出来。

for i indata1 data2 data3 data4 data5 data6

htseq-count-s no -t miRNA -i ID -o $i.hc.sam $i.ht.sammiRNA_reference/miRNA.gff3 | tee$i.count &

ls$i.count>> count.list

done

第六，然后就是重头戏，差异表达的miRNA，这里分为有重复的处理和没有重复的处理两种，对于没有重复的用DEGseq处理，有重复的用DESeq处理。

没有重复的用DEGseq：

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)