查看: 354|回复: 9

[流程分析] 玉米RNA-seq测序数据差异基因分析

  [复制链接]
  • TA的每日心情
    好棒
    2018.11.23 10:02
  • 签到天数: 112 天

    连续签到: 1 天

    [LV.6]常住居民II

    钵水母

    Rank: 3Rank: 3

    主题
    2
    奥币
    1124
    积分
    80
    注册时间
    2016.9.13
    在线时间
    38 小时

    发表于 2018.10.29 17:19:17 | 显示全部楼层 |阅读模式
    软件及参考基因组:
    BWA, Samtools, Hisat2, HTseq, gffcompare, Stringtie, BallgownR
    B73-V3, B73-V4(Enzampl database)
    流程:protocol首先从原始RAN-seq数据入手,先经过质控fastqc,之后检测rRNA占比,去除杂的reads之后进行数据处理;使用HISAT2将读段匹配到参考基因组上,可提供注释文件;StringTie进行转录本组装,估算每个基因及isoform的表达水平;所有的转录本merge的数据再一次被呈递给StringTie,重新估算转录本的丰度,提供转录本reads数量的数据给下一步的ballgownBallgown从上一步获得所有转录本及其丰度,根据实验条件进行分类统计。
    具体操作:
    1. 质控检测
    fastqc *.gz
    Fail选项:
    1Per Base Sequence Content(对所有reads的每一个位置,统计ATCG四种碱基正常情况的分布):当任一位置的A/T比例与G/C比例相差超过20%,报"FAIL"
    2Per Sequence GC Content:偏离理论分布的reads超过30%时,报"FAIL"
    3Duplicate Sequences统计序列完全一样的reads的频率。测序深度越高,越容易产生一定程度的duplication,这是正常的现象,但如果duplication的程度很高,就提示我们可能有bias的存在(如建库过程中的PCR duplication):当非uniquereads占总数的比例大于50%时,报"FAIL"
    2. rRNA检测
    NCBI网站下载Zea mays Eukaryotic 5S ribosomal RNAZea mays Eukaryotic 18S ribosomal RNAZea mays Eukaryotic 28S ribosomal RNA
    # bwa index rRNA-maize.fa
    # bwa mem -t 8 ../ref/rRNA-maize.fa NCSU_Wang_120430HiSeqRun_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   
    3. 玉米参考基因组及注释下载
    # 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
    从注释文件里面抽取出剪切位点和外显子信息,建立Index:
    # extract_splice_sites.py Zea_mays.AGPv4.38.gtf > maize-v4.ss
    #extract_exons.py Zea_mays.AGPv4.38.gtf > maize-v4.exon
    # hisat2-build --ss maize-v4.ss --exon maize-v4.exon Zea_mays.AGPv4.dna.toplevel.fa maize-v4_tran 速度慢,限速
    4.hisat2比对
    #hisat2 -p 16 -x ./maize/Zea_mays_tran -U ./48/NCSU_Wang_120430HiSeqRun_B73-mt-48h.fastq.gz -S ./result/B73.sam
    ###多组数据用逗号隔开,双末端数据用-1-2,单末端数据用-U,可以两种参数同时使用,reads长度可以不一致
    # python -m HTSeq.scripts.count 统计reads number
    5. 组装转录本并定量表达基因
    # stringtie B73.bam -p 8 -G
    多个样品:
    for i in *.sam;
    do
    i=${i%.sam*}; nohup samtools sort -@ 8 -o ${i}.bam ${i}.sam &
    done
    6.差异基因分析
    R
    >install.packages("devtools",repos="http://cran.us.r-project.org")
    >source("http://www.bioconductor.org/biocLite.R")
    >biocLite("ballgown", "RSkittleBrewer", "genefilter","dplyr","devtools")
    >library(ballgown)
    >library(RSkittleBrewer)
    >library(genefilter)
    >library(dplyr)
    >library(devtools)
    #  vim phenodata.csv
    ids, phenotype
    sample-s1B73mt-22,mt
    sample-s1B73WT-22,WT
    sample-s2B73mt-22,mt
    sample-s2B73WT-22,WT
    sampleB73mt-22,mt
    sampleB73WT-22,WT
    # vim mergelistb73-22.txt
    B73WT-22.gtf
    B73mt-22.gtf
    B73-WT-22.gtf
    B73-mt-22.gtf
    B73-WT-22.gtf
    B73-mt-22.gtf
    # stringtie --merge -p 8 –G Zea_mays.AGPv3.31.gtf -o b73_22.gtf mergelist_v3.txt
      stringtie -e -B -p 8 -G b73_22.gtf –o sampleB73WT-22.gtf sampleB73WT-22.bam
      stringtie -e -B -p 8 -G b73_22.gtf –o sampleB73mt-22.gtf sampleB73mt-22.bam
    stringtie -e -B -p 8 -G b73_22.gtf –o sample-s1B73WT-22.gtf sample-s1B73WT-22.bam
    stringtie -e -B -p 8 -G b73_22.gtf –o sample-s1B73mt-22.gtf sample-s1B73mt-22.bam
    stringtie -e -B -p 8 -G b73_22.gtf –o sample-s2B73WT-22.gtf sample-s2B73WT-22.bam
    stringtie -e -B -p 8 -G b73_22.gtf –o sample-s2B73mt-22.gtf sample-s2B73mt-22.bam
    7.差异基因作图
    7.1 导入数据
    >pheno_data = read.csv("phenodata.csv")
    Note: The sample names you have stored in your phenotype file do not match the file names of the samples you ran with StringTie!!!
    > all(pheno_data$ids == list.files("ballgown"))
    检查ballgownlist顺序:
    list.files(“ballgown”)
    bg = ballgown(dataDir = "ballgown",samplePattern = "sample",pData=pheno_data)
    bg = ballgown(dataDir = "ballgown", samplePattern = "sample", meas = "all")
    7.2滤掉低丰度的基因
    >bg_filt = subset(bg,"rowVars(texpr(bg))>1",genomesubset=TRUE)
    7.3确认组间有差异的转录本
    >result_transcripts = stattest(bg_filt,feature="transcript",covariate="phenotype",getFC=TRUE,meas="FPKM")
    > result_transcripts
    7.4确认组间有差异的基因
    >result_gene = stattest(bg_filt,feature="gene",covariate="phenotype",getFC=TRUE,meas="FPKM")
    > result_gene
    7.5 对结果 result_transcripts增加基因名
    >result_transcripts =data.frame(geneNames=ballgown::geneNames(bg_filt),geneIDs=ballgown::geneIDs(bg_filt), result_transcripts)
    7.6按照P值排序(从小到大)
    > result_transcripts=arrange(result_transcripts,pval)
    > result_gene=arrange(result_gene,pval)
    7.7把结果写入csv文件
    > write.csv(result_transcripts, "transcript_results.csv",row.names=FALSE)
    > write.csv(result_gene, "gene_results.csv",row.names=FALSE)
    7.8筛选出q值小于0.05transcriptsgenes
    >subset(result_transcripts,result_transcripts$qval<0.05)
    >subset(result_gene,result_gene$qval<0.05)
    > write.csv(bg, "bg",row.names=FALSE)
    8. 数据可视化之颜色设定
    # 赋予调色板五个指定颜色 > tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow') > palette(tropical)
    # 当然rainbow()函数也可以完成这个任务 > palette(rainbow(5))
    10.1FPKM为参考值作图,以time作为区分条件
    # 提取FPKM
    > fpkm = texpr(bg,meas="FPKM")
    #方便作图将其log转换,+1是为了避免出现log2(0)的情况 > fpkm = log2(fpkm+1)
    # 作图 > boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
    图1
    10.2就单个转录本的查看在样品中的分布
    > ballgown::transcriptNames(bg)[12]
    12 "NR_027232"
    > ballgown::geneNames(bg)[12]
    12 "LINC00685" #
    绘制箱线图
    > plot(fpkm[12,] ~ pheno_data$phenotype, border=c(1,2),
    + main=paste(ballgown::geneNames(bg)[12],' : ',
    + ballgown::transcriptNames(bg)[12]),pch=19, xlab="phenotype",
    + ylab='log2(FPKM+1)')    ## (左图)
    > points(fpkm[12,] ~ jitter(as.numeric(pheno_data$phenotype)),
    + col=as.numeric(pheno_data$phenotype))  ## (右图)
    10.3查看某一基因位置上所有的转录本
    >plotTranscripts(ballgown::geneIDs(bg)[442], bg, main=c('Gene MSTRG.168 in sample sampleB73WT-22'), sample=c('sampleB73WT-22'))
    有空格就出错 图2
    # plotTranscripts函数可以根据指定基因的id画出在特定区段的转录本
    #可以通过sample函数指定看在某个样本中的表达情况,这里选用id=442, sample= sampleB73WT-22
    > plotMeans('MSTRG.168', bg_filt,groupvar="phenotype",legend=FALSE)
    图3
















    本帖子中包含更多资源

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x

    评分

    参与人数 2奥币 +31 贡献 +2 收起 理由
    zhouyulu + 16 + 2 Well done!
    基迪奥-李泽标 + 15 楼主V5!

    查看全部评分

    回复

    使用道具 举报

  • TA的每日心情
    吃饭
    7 小时前
  • 签到天数: 335 天

    连续签到: 11 天

    [LV.8]以坛为家I

    版主

    Rank: 10Rank: 10Rank: 10

    主题
    60
    奥币
    6320
    积分
    2657
    注册时间
    2017.9.21
    在线时间
    231 小时

    突出贡献优秀版主


    发表于 2018.10.29 19:28:56 | 显示全部楼层
    nice!请问楼主是中农大玉米中心的吗?
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    忙~
    1 小时前
  • 签到天数: 296 天

    连续签到: 7 天

    [LV.8]以坛为家I

    帝王蝶

    Rank: 4

    主题
    0
    奥币
    1057
    积分
    310
    注册时间
    2016.11.16
    在线时间
    125 小时

    发表于 2018.10.30 08:40:25 | 显示全部楼层
    非常棒的流程,感谢分享
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    yes!
    2018.10.30 09:22
  • 签到天数: 2 天

    连续签到: 1 天

    [LV.1]初来乍到

    草履虫

    Rank: 2

    主题
    0
    奥币
    368
    积分
    8
    注册时间
    2016.5.11
    在线时间
    1 小时

    发表于 2018.10.30 09:28:37 | 显示全部楼层
    非常棒的分享!
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    yes!
    昨天 13:12
  • 签到天数: 5 天

    连续签到: 1 天

    [LV.2]偶尔看看I

    草履虫

    Rank: 2

    主题
    0
    奥币
    119
    积分
    46
    注册时间
    2017.7.7
    在线时间
    18 小时

    发表于 2018.10.30 15:41:42 | 显示全部楼层
    厉害!!!
    回复

    使用道具 举报

  • TA的每日心情
    no
    3 天前
  • 签到天数: 120 天

    连续签到: 6 天

    [LV.7]常住居民III

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    467
    积分
    78
    注册时间
    2017.3.29
    在线时间
    40 小时

    发表于 2018.10.31 08:11:46 | 显示全部楼层
    厉害
    回复

    使用道具 举报

  • TA的每日心情

    4 天前
  • 签到天数: 235 天

    连续签到: 1 天

    [LV.7]常住居民III

    帝王蝶

    Rank: 4

    主题
    2
    奥币
    560
    积分
    384
    注册时间
    2016.3.1
    在线时间
    89 小时

    发表于 2018.11.1 21:09:41 | 显示全部楼层
    回复

    使用道具 举报

  • TA的每日心情

    2018.11.30 19:29
  • 签到天数: 114 天

    连续签到: 1 天

    [LV.6]常住居民II

    帝王蝶

    Rank: 4

    主题
    1
    奥币
    504
    积分
    282
    注册时间
    2016.9.2
    在线时间
    50 小时

    发表于 2018.11.5 18:40:53 | 显示全部楼层
    太棒了,感谢慷慨的楼主!
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    害羞
    6 小时前
  • 签到天数: 390 天

    连续签到: 2 天

    [LV.9]以坛为家II

    中华鲟

    Rank: 5Rank: 5

    主题
    15
    奥币
    1583
    积分
    819
    注册时间
    2016.1.8
    在线时间
    198 小时

    发表于 2018.11.28 10:39:49 | 显示全部楼层
    没用到gffcompare
    回复 支持 反对

    使用道具 举报

  • TA的每日心情

    3 天前
  • 签到天数: 81 天

    连续签到: 3 天

    [LV.6]常住居民II

    钵水母

    Rank: 3Rank: 3

    主题
    0
    奥币
    264
    积分
    114
    注册时间
    2018.8.28
    在线时间
    26 小时

    发表于 2018.11.28 11:00:39 | 显示全部楼层
    收藏
    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    快速回复 返回顶部 返回列表