b73 玉米RNA-seq测序数据差异基因分析

今天欢鹰给大家分享了一个很棒的玉米转录组流量分析。原作者是cxge，最早发表在omicshare论坛。看完原文，可以跳转到这篇文章的帖子~

软件和参考基因组

BWA，萨姆托勒，Hisat2，HTseq，gffcompare，Stringtie，ball Loma，R

B73-V3，B73-V4(Enzampl数据库)

流动

协议从原始RAN-seq数据开始，先经过质量控制fastqc，然后检测rRNA的比例，再去除杂读进行数据处理；使用HISAT2将阅读片段与参考基因组匹配，可以提供注释文件；StringTie组装转录物以估计每个基因和同工型的表达水平；所有转录合并的数据再次呈现给StringTie，转录的丰度被重新估计，并且转录读取的数量被提供给下一个ball loat；；ball loat从上一步获得了所有的转录本及其丰度，并根据实验条件进行了分类。

具体操作

一个

质量控制测试

fastqc *。地面零点

失败选项:

(1)每个碱基序列内容(对于所有读数的每个位置，计算ATCG四个碱基的正态分布):当任意位置的A/T比和G/C比之差超过20%时，报告“失败”。

(2)每序列气相色谱含量:当偏离理论分布的读数超过30%时，将报告“失败”。

(3)重复序列计算具有相同序列的读取频率。测序深度越高，越容易产生一定程度的重复，这是正常现象。但如果重复程度很高，则说明我们可能有偏差(比如数据库构建过程中的PCR重复):当非唯一性读取占总数的比例超过50%时，就会被报为“FAIL”

2

RRNA探测

从NCBI网站下载玉米真核5S核糖体RNA玉米真核18S核糖体RNA玉米真核28S核糖体RNA

# bwa指数RRna-玉米. fa

# bwa mem -t 8../ref/RRNA-玉米. fa _王_ 120430 his qrun _ B73-WT-22.fastq.gz | samtools view-bS->；B73-WT.bam

# samtools sort B73-WT . bam-o B73-WT . sort . bam

# samtools flagstat B73-WT.bam

三

玉米参考基因组和注释下载

# wget ftp://ftp.ensemblgenomes.org/pub/release-38/plants/fasta/ Zea _ mays/DNA/Zea _ mays。AGPv4.dna.toplevel.fa.gz

# wget ftp://ftp.ensemblgenomes.org/pub/release-31/plants/fasta/ Zea _ mays/DNA/Zea _ mays。AGPv3.31.dna.genome.fa.gz

# wget ftp://ftp.ensemblgenomes.org/pub/release-38/plants/gtf/ Zea _ mays/Zea _ mays。AGPv4.38.gtf.gz

# wget ftp://ftp.ensemblgenomes.org/pub/release-31/plants/gtf/ Zea _ mays/Zea _ mays。AGPv3.31.gtf.gz

从注释文件中提取切割位点和外显子信息，并建立索引:

# extract _ splice _ sites . py Zea _ mays。AGPv4.38.gtf >。玉米-v4.ss

#extract_exons.py Zea_mays。AGPv4.38.gtf >。玉米-v4 .外显子

# his T2-build-ss玉米-v4 . ss-外显子玉米-v4 .外显子Zea _ mays . agpv 4 . DNA . toplevel . famaise-v4 _ tran速度慢，速度受限

四

Hisat2对齐

# his T2-p 16-x ./玉蜀黍/Zea _ mays _ tran-U ./48/_ Wang _ 120430 Hiseqrun _ B73-mt-48h . fastq . gz-S ./result/B73 . Sam

# # #多组数据用逗号分隔，双端数据为-1和-2，单端数据为-U，可以同时使用，读取长度可以不一致

# python-m HTSEq . s . count statistics读取数量

五

组装转录物并定量表达基因

# stringtie B73.bam -p 8 -G

多个样本:

因为我在。萨姆。

做

i=${i%。Sam * }；nohup samtools sort -@ 8 -o ${i}。bam ${i}。sam & amp

行

六

差异基因分析

稀有

>。install.packages("devtools "，repo = " http://cran . us . r-project . org "/>

8.2检查样本中单个转录本的分布

>。ball loat::tran names(BG)[12]

12 "NR_027232 "

>。* GeneName(BG)[12]

12 "LINC00685 "

#画一个方框图

>。图(fpkm[12，]~ pheno _ data $显型，border=c(1，2)、

+main = paste(ball袍::geneNames(bg)[12]，' : '，

+ball长袍::tranNames(bg)[12])，pch=19，xlab= "显型"，

+ ylab='log2(FPKM+1)')

# #(左图)

>。点数(fpkm[12，] ~抖动(as . numeric(pheno _ data $显型))，

+col = as . numeric(pheno _ data $显型))

# #(右)

8.3检查某个基因位置的所有转录本

>。plotTrans(ball袍::geneIDs(bg)[442]，bg，main = c(' Gene MSTRG . 168 in sample b73wt-22 ')，sample=c('sampleB73WT-22 ')

空网格出错图2

# plotTrans函数可以根据指定基因的id绘制特定片段中的转录本

#您可以通过sample函数指定样本中的表达式，其中选择了id = 442和sample = sampleb73wt-22

>。plotMeans('MSTRG.168 '，bg_filt，groupvar= "显型"，图例=FALSE)

图3

今天就到这里~。更多圣信干货，可以登录omisshare论坛(http://www.omicshare.com/forum/)阅读原文或直接跳转到omisshare论坛。