• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)跟着豆豆一起探索可变剪切分析

      (未测试)跟着豆豆一起探索可变剪切分析

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2020年2月23日
      • 评论 0评论

      介绍

      概念介绍

      可变剪接(Alternative splicing)是基因表达的一种调节方式。真核细胞的基因序列中,包含了内含子与外显子,两者交互穿插。一般来说,内含子在基因转录成mRNA前体后会被RNA剪接体移除,剩下的外显子才是能够存在于成熟mRNA(之后再进一步转译成蛋白质)的片段。一条未经剪接的RNA,含有的多种外显子被剪成的不同组合,可转译出不同的蛋白质。这样保证了:同一基因中的外显子以不同形式进行组合,使一个基因在不同时间、不同环境中能够制造出不同的蛋白质,增加生理状况下系统的复杂性或适应性。

      其他解释:

      1:有些基因的一个mRNA前体通过不同的剪接方式(选择不同的剪接位点)产生不同的mRNA剪接异构体(isoform),这一过程称为可变剪接(或选择性剪接,alternative splicing)

      2:转录组一般是指从细胞或组织的基因组所转录出来的RNA的总和,包括编码蛋白质的mRNA和各种非编码RNA(rRNA,tRNA,snRNA,snoRNA,lncRNA,microRNA等)。真核生物的基因结构是不连续的。基因组最初的转录产物其实并不是成熟的mRNA分子,而是它的前体pre-mRNA,那么怎么变成成熟的mRNA呢?

      就需要从pre-mRNA中将非编码蛋白质的内含子(intron)切除,然后拼接剩下的编码蛋白质的外显子(exon)。但实际上,在这个过程中,有多种多样的前切和拼接方式,从而产生不同的剪切异构体,即可变剪切事件。

      施一公在2017未来科学大奖颁奖典礼发表主题演讲:剪接过程如此之复杂,你需要七个不同的剪接体,准确的识别内含子,把外显子拼在一起,这个过程很容易出错。目前人类35%的遗传疾病或由RNA选择性剪接导致,比如视网膜色素变性、脊髓性肌肉萎缩症等。目前已有针对剪接体的药物被美国食品药品检验局批准上市,价格昂贵但非常有效
      http://tech.sina.com.cn/2017-10-29/doc-ifynhhay8118096.shtml

      人基因组中,大约95%的多外显子基因存在可变剪接。单独一个基因通过可变剪接产生十几种剪接异构体的现象很常见。最突出的例子是果蝇的Dscam基因,其潜在的可变剪接类型有38016个。

      与疾病相关的几种典型剪接变异:

      1) 点突变–地中海贫血

      2) 与疾病相关的短重复元件突变–Alport综合症

      3) 单核苷酸多态性–阿兹海默症

      4) 可变剪接产生蛋白质亚型的比例改变–神经退化性疾病

      5) 剪接因子的缺失–脊髓型肌萎缩

      研究进展:https://www.nature.com/subjects/alternative-splicing

      一篇介绍AS机制的文献:Mechanism of alternative splicing and its regulation

      类型介绍

      可变剪接是调节基因表达和产生蛋白质组多样性的重要机制,是导致真核生物基因和蛋白质数量较大差异的重要原因。在生物体内,主要存在7种可变剪接类型:

      • ES(Exon skipping):外显子跳跃。外显子在前体mRNA剪接形成成熟mRNA过程中被跳过,最终没有出现在某些成熟mRNA上【发生跳跃的外显子和其两侧的内含子都被剪切掉;上游和下游的外显子被直接连着一起保留在剪切后的产物中】
      • RI(Retained intron):内含子保留。前体mRNA在剪接形成成熟mRNA的过程中,部分内含子被保留下来【某一段核苷酸序列在一个剪切体中是外显子的一部分,而在与之对照的剪切体中却是内含子而被剪切掉】
      • AD(Alternate Donor site)或A5SS(Alternative 5′ splice site):5’端可变剪接。前体mRNA在剪接形成成熟mRNA的过程中,5’端边界发生不同方式的剪接,导致5‘端外显子有所延长
      • AA(Alternate acceptor site)或A3SS(Alternative 3′ splice site):3’端可变剪接。前体mRNA在剪接形成成熟mRNA的过程中,3’端边界发生不同方式的剪接,导致3‘端外显子有所延长
      • AT(Alternate terminator)或Alternative last exon:第一个外显子发生改变
      • AP(Alternate promoter)或Alternative first exon:最后一个外显子发生改变
      • ME(Mutually exclusive exon):外显子选择性跳跃。形成的两种不同的转录本中,两转录本之间相同的外显子称为constitutive exon, 不同的外显子称为inclusive exon,inclusive exon不能同时存在与同一转录本中, 只能分别存在于不同转录本中

      分析工具

      1=> rMATS

      介绍

      网站:http://rnaseq-mats.sourceforge.net/

      目前最新版是4.0.2

      之前叫MATS(Multivariate Analysis of Transcript Splicing),后来为了有生物重复的RNA-Seq又升级到rMATS,相比之前的3.2.5版本,计算速度提升,并可支持多线程执行

      它能识别5种AS事件

      安装

      默认官网安装十分复杂:需要检测各种python,还需要sudo apt-get等等root权限。但十有八九最后还是会遇到一个问题:error while loading shared libraries: libgsl.so.0: cannot open shared object file: No such file or directory。此时依然是需要root安装一个库

      • 使用:sudo apt-get install libgsl0-dev
      • 解决方案依然是需要root:https://www.biostars.org/p/336864/

      直接conda安装:

      conda install -c bioconda rmats rmats2sashimiplot star

      安装好检测:

      rmats.py --help

      用法说明如下:

      使用示例

      下载官网的示例数据:

      mkdir ~/biosoft/rmats_v4.0.2 && cd ~/biosoft/rmats_v4.0.2
      # -O 指定以不同文件名保存
      wget https://sourceforge.net/projects/rnaseq-mats/files/MATS/testData.tgz/download -O testData.tgz
      # 然后解压
      tar xvzf testData.tgz
      
      4.5M Aug 15  2016 231ESRP.25K.rep-1.R1.fastq
      4.5M Aug 15  2016 231ESRP.25K.rep-1.R2.fastq
      2.6M Apr 26  2018 231ESRP.25K.rep-1.bam
      4.6M Aug 15  2016 231ESRP.25K.rep-2.R1.fastq
      4.6M Aug 15  2016 231ESRP.25K.rep-2.R2.fastq
      2.6M Apr 26  2018 231ESRP.25K.rep-2.bam
      4.5M Aug 15  2016 231EV.25K.rep-1.R1.fastq
      4.5M Aug 15  2016 231EV.25K.rep-1.R2.fastq
      2.6M Apr 26  2018 231EV.25K.rep-1.bam
      4.6M Aug 15  2016 231EV.25K.rep-2.R1.fastq
      4.6M Aug 15  2016 231EV.25K.rep-2.R2.fastq
      2.6M Apr 26  2018 231EV.25K.rep-2.bam
        44 Apr 26  2018 b1.txt
        40 Apr 26  2018 b2.txt
       108 Apr 26  2018 s1.txt
       100 Apr 26  2018 s2.txt

      软件支持两种格式输入,文件内以逗号分隔重复样本名:

      • 使用fastq测序数据:这样的话需要用--s1、--s2分别指定两个样本的fq文件名,这样方便了自己,不需要自己运行STAR比对,直接用fq文件出结果。但需要自己下载安装STAR和比对的索引文件

      Obtain STAR genome index for genome by either of the following two ways:

      • Download pre-built STAR indexes if using Human (hg38, hg19) or Mouse (mm10). 【hg38、hg19、mm10有现成的可以下载:http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STARgenomes/】
      • Build your own STAR index following STAR manual from genome fasta sequence【自己构建的话:
        https://github.com/alexdobin/STAR/tree/master/doc】
      • 使用比对后的bam:它支持其他比对软件(如tophat、hisat2等)作为输入,从而减少rMATS的运行时间
      如果直接使用fastq测序数据
      # 测试数据中已经帮我们准备好了名称
      cat s1.txt
      231ESRP.25K.rep-1.R1.fastq:231ESRP.25K.rep-1.R2.fastq,231ESRP.25K.rep-2.R1.fastq:231ESRP.25K.rep-2.R2.fastq
      
      cat s2.txt
      231EV.25K.rep-1.R1.fastq:231EV.25K.rep-1.R2.fastq,231EV.25K.rep-2.R1.fastq:231EV.25K.rep-2.R2.fastq
      

      我们自己需要准备好GTF文件和STAR index,然后:

      # 打开下面?链接
      # GTF和STAR for hg19:
      http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STARgenomes/ENSEMBL/homo_sapiens/ENSEMBL.homo_sapiens.release-75/
      
      # GTF和STAR for hg38:
      http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STARgenomes/ENSEMBL/homo_sapiens/ENSEMBL.homo_sapiens.release-83/
      
      # 如果运行STAR
      STAR --readFilesIn test_1.clean.fq test_2.clean.fq \
      --chimSegmentMin 2 \
      --outFilterMismatchNmax 3 \
      --alignEndsType EndToEnd \
      --runThreadN 5 \
      --outSAMstrandField intronMotif \
      --outSAMtype BAM SortedByCoordinate \
      --alignSJDBoverhangMin 6 \
      --alignIntronMax 299999 \
      --genomeDir /YOUR_INDEX_DIR/ \
      --sjdbGTFfile /YOUR_GTF_DIR/ \
      --outFileNamePrefix test 
      
      # 然后运行rmats
      rmats.py --s1 s1.txt --s2 s2.txt --gtf gtf/Homo_sapiens.Ensembl.GRCh37.72.gtf --bi ~/STARindex/hg19 --od out_test -t paired --nthread 6 --readLength 50 --tophatAnchor 8 --cstat 0.0001 --tstat 6
      如果直接使用bam测序数据
      # 测试数据中已经帮我们准备好了名称
      cat b1.txt
      231ESRP.25K.rep-1.bam,231ESRP.25K.rep-2.bam
      
      cat b2.txt
      231EV.25K.rep-1.bam,231EV.25K.rep-2.bam

      只需要提供GTF即可:

      rmats.py --b1 b1.txt --b2 b2.txt --gtf ~/reference/gtf/ensembl/hg19/Homo_sapiens.GRCh37.75.gtf --od bam_test -t paired --readLength 50 --cstat 0.0001 --libType fr-unstranded
      看一些参数的意思
      -t readType 双端测序为paired,单端测序为single
      --readLength 测序reads的长度(这里测试数据的长度是51)
      --gtf 指定GTF文件
      --od 指定输出目录
      --nthread 设置线程数
      --cstat 设置一小数 (the cutoff splicing difference,default is 0.0001 for 0.01% difference. Valid: 0 ≤ cutoff < 1)【它采用的是Likelihood-ratio test来比较两组样本可变剪切是否有差异】
      -libType 默认是非链特异性(fr-unstranded);如果使用特异性文库,可以指定fr-firststrand或者fr-secondstrand
      结果文件

      结果的解读看:https://cloud.tencent.com/developer/article/1366294

      可视化

      使用rmats2sashimiplot

      需要注意:

      • 如果直接使用rMATS从fastq输入文件到结果输出文件,那么rMATS是不保留中间的bam文件的
      • rmats2sashimiplot可视化则需要bam文件作为输入,因此最好自己先进行一轮比对,然后进行rMATS,最后再plot
      • 使用的bam文件要先用samtools对bam文件建索引,否则rmats2sashimiplot会自己再做一遍,而且很慢
      • 确保rMATS的结果文件中染色体和bam文件的染色体名称一致(要么都是chr1,要么都是直接数字1)
      • 出图结果默认是PDF
      # 以单端为例
      rmats2sashimiplot --b1 A1.bam,A2.bam,A3.bam \
      --b2 B1.bam,B2.bam,B3.bam \
      -t SE \
      -e ./SE.MATS.JC_top20.txt \
      --l1 A --l2 B \
      --exon_s 1 --intron_s 1 \
      -o SE_plot

      2 => CircSplice

      介绍
      http://gb.whu.edu.cn/CircSplice/

      circRNA作为一个新兴的非编码RNA,也同样通过mRNA前体的可变剪切而来。作者建议:测序数据文库最好采用去rRNA及线性RNA。

      目前常用的可变剪切事件预测软件有 CIRI-full (更适合长片段,>250或300bp),以及短片段的CIRI-AS (50~150bp),但是并没有提供样本间比较的功能。

      它是一款基于perl编写的流程,能够通过预测反向剪切事件识别circRNA,支持 GT-AG和CT-AC两种剪切位点。提供四种circ-AS事件类型结果:外显子跳跃(SE),内含子保留(RI),5’端可变剪接(A5SS)和3’端可变剪接(A3SS)

      1 质控及过滤
      fastqc Sample.R1.fq.gz Sample.R2.fq.gz
      trim_galore -q 20 --phred33 --stringency 3 --length 20 -e 0.1 --paired Sample.R1.fq Sample.R2.fq --gzip -o Sample
      fastqc Sample.R1_trimmed.fq.gz Sample.R2_trimmed.fq.gz
      2 比对
      # 建索引
      STAR --runThreadN NumberOfThreads --runMode genomeGenerate --genomeDir /path/to/genomeDir --genomeFastaFiles /path/to/genome/fasta1 /path/to/genome/fasta2 --sjdbGTFfile /path/to/annotations.gtf
      # 比对
      STAR --genomeDir /path/to/genomeDir --readFilesIn Sample.R1_trimmed.fq.gz Sample.R2_trimmed.fq.gz --readFilesCommand zcat --runThreadN 10 --chimSegmentMin 20 --chimScoreMin 1 --alignIntronMax 100000 --outFilterMismatchNmax 4 --alignTranscriptsPerReadNmax 100000 --outFilterMultimapNmax 2 --outFileNamePrefix Sample
      3 鉴定
      # For paired-end library,
      CircSplice.pl Chimeric.out.sam hg38.genome.fa bed_refFlat_hg38.txt
      
      # single-end library,
      CircSplice-single.pl Chimeric.out.sam hg38.genome.fa bed_refFlat_hg38.txt

      将获得两个结果文件,Chimeric.out.sam.result.as(circ-As时间)和Chimeric.out.sam.result.circ(circRNA鉴定)

      4 合并结果
      CircSplice-merge.pl dir-as dir-circ

      CircSplice-merge merges results according to the genomic coordinates with 2bp mismatch toleration.

      得到结果

      详细的结果见:http://gb.whu.edu.cn/CircSplice/userguide.html

      • AS_result 中每一行代表样本中的一个AS事件
      • Circ_result 中每一行代表一条circRNA
      其他物种

      CircSplice.pl 也支持其他物种或者基因组的版本,需要自行构建,文本格式如下

      Chromosome
      Start coordinates -2
      End coordinates +2
      Transcript ID
      Number of exons
      Strand
      Gene symbol
      Transcript ID
      Chromsome
      Strand
      Start coordinates
      end coordinates
      Start coordinates of CDS
      End coordinates of CDS
      Number of exons
      Start coordinates of each exon
      End coordinates of each exon
      Gene type (lncRNA or mRNA)

      再通过 reftobed.pl refFlat.txt转换bed-refFlat格式

      3 => CASH

      介绍
      http://www.novelbio.com/cash/index.html

      这是一款国产软件,全称:Comprehensive alternative splicing hunting

      It aims to self-construct AS (alternative splicing) sites and detect differential AS events between samples of RNA-Seq data.
      引用:
      Zhou X, Wu W, Li H, Cheng Y, Wei N, Zong J, Feng X, Xie Z, Chen D, and Manley JL et al. 2014. Transcriptome analysis of alternative splicing events regulated by SRSF10 reveals position-dependent splicing modulation. Nucleic Acids Res 42: 4019-4030.

      包含了主要的两步:  SpliceCons (Splice site Construction) and SpliceDiff(differential AS detection)

      • SpliceCons increases the recognition of AS events considerably and subsequently
      • SpliceDiff uses two combined statistical methods to improve the detection of differential AS events

      需要提供的文件:比对产生的bam文件(必须经过sort + index)、GTF基因注释文件(另外,如果想找到剪切位点附近的序列,还需要上传基因组序列信息)

      使用说明在:http://www.novelbio.com/cash/manual_index.html

      下载地址在:http://www.novelbio.com/cash/cash/cash_v2.2.0.zip

      参考数据在:http://www.novelbio.com/cash/cash/Examples.zip

      下载好检查软件
      unzip cash_v2.2.0.zip
      cd cash_v2.2.0
      # 需要提前用conda install -c anaconda openjdk -y
      java -jar cash.jar --help
      使用示例
      java -jar -Xmx10g cash.jar \
      --Case:Mutation file1.bam,file2.bam \
      --Control:WildType file3.bam,file4.bam \
      --GTF genes.gtf --Output outFilePrefix

      需要JAVA环境:requires jre1.8 or later, and at least 8GB memory for 2vs2 human samples

      关于结果

      4 => 【R包】LeafCutter v0.2.8

      文章在:Annotation-free quantification of RNA splicing using LeafCutter
      特点:可省略GTF/GFF文件;可以检测splicing quantitative trait loci (sQTLs) 官网在:http://davidaknowles.github.io/leafcutter/
      参考:http://www.bio-info-trainee.com/2949.html

      安装可以用

      conda install -c davidaknowles r-leafcutter
      # 或者直接在Rstudio
      devtools::install_github(“davidaknowles/leafcutter/leafcutter”)

      官方提供的一些脚本和测试数据在
      git clone https://github.com/davidaknowles/leafcutter
      cd leafcutter
      ## 其中有一个bam2junc.sh 一会要用到
      首先:bam2junc
      # $bam是bam名称; $sample.junc是生成的结果
      bam2junc.sh  $bam $sample.junc
      # 多个bam可以批量操作
      # 会得到
      chr11    109735445   109827551   .   19  +
      chr18    48458730    48465939    .   8   -
      chr12    82751048    82752457    .   12  -
      chr15    51018323    51018517    .   14  -
      然后:intron clustering

      需要python 2.7

      ls *.junc >test_juncfiles.txt
      
      python leafcutter_cluster.py -j test_juncfiles.txt -m 50 -o test -l 500000
      # 结果的**_perind_numers.counts.gz文件中,每一行是一个内含子,每一列是一个样本的表达值(然后利用这些数值进行可变剪切),例如:
      chr8:145651155:145651305:clu_6538 21 14 19 8 9 0 13 33 0 0 4 0 5 8 12 0 12 34 15 0 0 10 11
      接着:制作分组文件进行差异分析

      分组文件就是两列:每个样本属于哪个组(有且只能有2个分组;如果分组很多可以多次处理多次运行),例如:

      test1.bam control
      test2.bam control
      test3.bam trt
      test4.bam trt

      有了分组文件再运行

      leafcutter_ds.R --num_threads 4 \
       --exon_file=~/leafcutter/data/gencode19_exons.txt.gz \
      test_perind_numers.counts.gz group_info.txt

      会得到:leafcutter_ds_cluster_significance.txt 【这个文件需要去理解】

      最后:可视化
      ds_plots.R -e  ~/leafcutter/data/gencode19_exons.txt.gz test_perind_numers.counts.gz   group_info.txt leafcutter_ds_cluster_significance.txt -f 0.05

      5 => 【R包】SGSeq

      Splice event prediction and quantification from RNA-seq data
      完全基于R语言处理
      官方说明:
      https://bioconductor.org/packages/release/bioc/vignettes/SGSeq/inst/doc/SGSeq.html

      参考:http://www.bio-info-trainee.com/2890.html

      安装
      if (!requireNamespace("BiocManager", quietly = TRUE))
          install.packages("BiocManager")
      
      BiocManager::install("SGSeq")
      需要的数据之bam

      输入文件是bam(可以用hisat2、star等先比对好),然后构建好索引

      当然,安装好这个包以后,也会有自带的测试数据

      需要的数据之注释文件

      通过convertToTxFeatures()函数把 GRanges对象转化成了一个TxFeatures对象

      # for hg19 (这里使用chr16)
      library(TxDb.Hsapiens.UCSC.hg19.knownGene)
      txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
      txdb <- keepSeqlevels(txdb, "chr16")
      seqlevelsStyle(txdb) <- "NCBI"
      
      # To work with the annotation in the SGSeq framework, transcript features are extracted from the TxDb object using function convertToTxFeatures()
      txf_ucsc <- convertToTxFeatures(txdb)
      txf_ucsc <- txf_ucsc[txf_ucsc %over% gr]
      head(txf_ucsc)
      type(txf_ucsc)
      head(txName(txf_ucsc))
      head(geneName(txf_ucsc))

      这样就能标记下面五种类型了:

      • J (splice junction)
      • I (internal exon)
      • F (first/5′-terminal exon)
      • L (last/5′-terminal exon)
      • U (unspliced transcript).

      再用convertToSGFeatures() 把 TxFeatures 转成 SGFeatures

      sgf_ucsc <- convertToSGFeatures(txf_ucsc)
      head(sgf_ucsc)

      SGFeatures对象包含了:

      • J (splice junction)
      • E (disjoint exon bin)
      • D (splice donor site)
      • A (splice acceptor site)
      开始分析
      sgfc_ucsc <- analyzeFeatures(si, features = txf_ucsc)
      sgfc_ucsc
      ## class: SGFeatureCounts 
      ## dim: 42 8 
      ## metadata(0):
      ## assays(2): counts FPKM
      ## rownames: NULL
      ## rowData names(0):
      ## colnames(8): N1 N2 ... T3 T4
      ## colData names(6): sample_name file_bam ... frag_length lib_size

      返回了一个SGFeatureCounts对象,在colData中存储了样本信息,rowRanges中存储了splice graph features,counts和FPKM存储了每个基因的每个外显子的表达量以及两个外显子中间夹着的内含子的表达情况。

      Counts for exons and splice junctions are based on structurally compatible fragments.
      FPKMs for splice graph features can be visualized with function plotFeatures

      # 检查结果
      colData(sgfc_ucsc)
      rowRanges(sgfc_ucsc)
      head(counts(sgfc_ucsc))
      head(FPKM(sgfc_ucsc))
      可视化其中一个基因
      df <- plotFeatures(sgfc_ucsc, geneID = 1)

      df <- plotFeatures(sgfc_ucsc, geneID = 1)

      复杂一点的可视化 based on de novo prediction
      sgfc_pred <- analyzeFeatures(si, which = gr)
      head(rowRanges(sgfc_pred))
      sgfc_pred <- annotate(sgfc_pred, txf_ucsc)
      head(rowRanges(sgfc_pred))
      df <- plotFeatures(sgfc_pred, geneID = 1, color_novel = "red") 

      图中顶部灰色部分是:exons and splice junctions predicted from RNA-seq data are consistent with transcripts in the UCSC knownGene table

      红色的是:an unannotated exon (E3, shown in red) was discovered from the RNA-seq data, which is expressed in three of the four normal colorectal samples (N2, N3, N4).

      其他展现形式
      # splice junctions
      plotFeatures(sgfc_pred, geneID = 1, include = "junctions")
      # exon bins
      plotFeatures(sgfc_pred, geneID = 1, include = "exons")
      # both exon bins and splice junctions
      plotFeatures(sgfc_pred, geneID = 1, include = "both")
      # 还能设置其他参数
      # ”toscale“ controls which parts of the gene model are drawn to scale.
      plotFeatures(sgfc_pred, geneID = 1, toscale = "gene")
      plotFeatures(sgfc_pred, geneID = 1, toscale = "exon")
      plotFeatures(sgfc_pred, geneID = 1, toscale = "none")
      # 还有,【per-base read coverages】 and 【splice junction counts】 can be visualized with function plotCoverage()
      par(mfrow = c(5, 1), mar = c(1, 3, 1, 1))
      plotSpliceGraph(rowRanges(sgfc_pred), geneID = 1, toscale = "none", color_novel = "red")
      for (j in 1:4) {
        plotCoverage(sgfc_pred[, j], geneID = 1, toscale = "none")
      }
      预测剪切异构体(splice variant)对蛋白编码的潜在影响

      依然是先得到SGVariants这个对象再进行处理

      # 鉴定Splice variant
      sgvc_pred <- analyzeVariants(sgfc_pred)
      sgvc_pred
      mcols(sgvc_pred)
      
      # 对Splice variant进行定量
      variantFreq(sgvc_pred)
      plotVariants(sgvc_pred, eventID = 1, color_novel = "red")
      
      # 最后解释Splice variant
      library(BSgenome.Hsapiens.UCSC.hg19)
      seqlevelsStyle(Hsapiens) <- "NCBI"
      vep <- predictVariantEffects(sgv_pred, txdb, Hsapiens)
      vep
      ##   variantID     txName geneName                        RNA_change
      ## 1         2 uc002fjv.3    79791         r.88_89ins88+1798_88+1884
      ## 2         2 uc002fjw.3    79791     r.412_413ins412+1798_412+1884
      ## 3         2 uc010vot.2    79791 r.-105_-104ins-105+1798_-105+1884
      ##   RNA_variant_type                              protein_change
      ## 1        insertion   p.K29_L30insRINPRVKSGRFVKILPDYEHMAYRDVYTC
      ## 2        insertion p.K137_L138insRINPRVKSGRFVKILPDYEHMAYRDVYTC
      ## 3        insertion                                         p.=
      ##   protein_variant_type
      ## 1   in-frame_insertion
      ## 2   in-frame_insertion
      ## 3            no_change

      结果的每行描述了:the effect of a particular splice variant on an annotated protein-coding transcript.

      最后

      SGSeq does not implement statistical tests for differential splice variant usage. However, existing software packages such as DEXSeq (Anders, Reyes, and Huber 2012) and limma (Ritchie et al. 2015) can be used for this purpose.

      6 => 【R包】DEXSeq 通过外显子差异分析推测可变剪切

      直接看:http://www.bio-info-trainee.com/bioconductor_China/software/DEXSeq.html
      帮助文档在:
      http://bioconductor.org/packages/release/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.html
      官方建议:it is important to align them to the genome, not to the transcriptome, and to use a splice-aware aligner (i.e., a short-read alignment tool that can deal with reads that span across introns) such as TopHat2 (Kim et al. 2013), GSNAP (Wu and Nacu 2010), or STAR (Dobin et al. 2013)

      首先安装包
      ## 首先安装包
      # first prepare BioManager on CRAN
      if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")
      if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
      if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
      
      # use BiocManager to install
      Biocductor_packages <- c('DEXSeq','pasilla')
      for (pkg in Biocductor_packages){
          if (! require(pkg,character.only=T) ) {
              BiocManager::install(pkg,ask = F,update = F)
              require(pkg,character.only=T) 
          }
      }
      
      suppressPackageStartupMessages(library("DEXSeq"))
      suppressPackageStartupMessages(library("pasilla"))
      准备3类输入文件=》1.counts矩阵

      可以利用DEXSeq包中的python脚本得到

      python /path/to/library/DEXSeq/python_scripts/dexseq_count.py Dmel.BDGP5.25.62.DEXSeq.chr.gff untreated1.sam untreated1fb.txt

      有了这个外显子count文件,再看一下:

      # 打开测试数据
      inDir <- system.file("extdata", package="pasilla")
      countFiles <- list.files(inDir, pattern="fb.txt$", full.names=TRUE)
      countFiles
      > ## 共七个样本,可以从文件名看到样本描述,即:每个基因的每个外显子的reads数量
      > head(read.table(countFiles[1]))
                     V1 V2
      1 FBgn0000003:001  0
      2 FBgn0000008:001  0
      3 FBgn0000008:002  0
      4 FBgn0000008:003  0
      5 FBgn0000008:004  1
      6 FBgn0000008:005  4
      准备3类输入文件=》2.gtf格式的基因组注释文件

      依然可以使用自带的python脚本制作:

      python /path/to/library/DEXSeq/python_scripts/dexseq_prepare_annotation.py Drosophila_melanogaster.BDGP5.72.gtf Dmel.BDGP5.25.62.DEXSeq.chr.gff

      看一下示例数据:

      gffFile <- list.files(inDir, pattern="gff$", full.names=TRUE) 
      gffFile
      准备3类输入文件=》3.实验设计表格
      sampleTable <- data.frame(row.names=c(paste("treated", 1:3, sep=""), paste("untreated", 1:4, sep="")),
                                condition=rep(c("knockdown", "control"), c(3, 4)))
      sampleTable
      ##            condition
      ## treated1   knockdown
      ## treated2   knockdown
      ## treated3   knockdown
      ## untreated1   control
      ## untreated2   control
      ## untreated3   control
      ## untreated4   control
      构造DEXSeqDataSet对象
      dxd <- DEXSeqDataSetFromHTSeq(
         countFiles,
         sampleData=sampleTable,
         design= ~sample + exon + condition:exon,
         flattenedfile=gffFile
         )
      
      > dxd
      class: DEXSeqDataSet 
      dim: 70463 14 
      metadata(1): version
      assays(1): counts
      rownames(70463): FBgn0000003:E001
        FBgn0000008:E001 ... FBgn0261575:E001
        FBgn0261575:E002
      rowData names(5): featureID groupID exonBaseMean
        exonBaseVar transcripts
      colnames: NULL
      colData names(3): sample condition exon
      
      > class(dxd)
      [1] "DEXSeqDataSet"
      attr(,"package")
      [1] "DEXSeq"
      差异分析
      dxr <- DEXSeq(dxd) # 耗时
      # 实际上做了以下几件事
      if(F){
          dxd <- estimateSizeFactors(dxd) 
          dxd <- estimateDispersions(dxd) 
          dxd <- testForDEU(dxd) 
          dxd <- estimateExonFoldChanges(dxd, fitExptoVar="condition")
          dxr1 <- DEXSeqResults(dxd)
      }
      指定基因进行可视化
      ## 最简单的图(红色部分是显著差异的外显子)
      plotDEXSeq( dxr2, "FBgn0010909", legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )
      
      ##将外显子表达量和转录本模型结合(putting differential exon usage results into the context of isoform regulation),这些转录本模型都是有注释的
      plotDEXSeq( dxr, "FBgn0000210", 
                  displayTranscripts=TRUE, legend=TRUE,
                  cex.axis=1.2, cex=1.3, lwd=2 )
      
      
      # 按表达量查看(counts are normalized by dividing them by the size factors)
      plotDEXSeq( dxr, "FBgn0000210", 
                  expression=FALSE, norCounts=TRUE,
                  legend=TRUE, cex.axis=1.2, 
                  cex=1.3, lwd=2 )

      最后还能汇总,得到所有的显著性差异的基因图片

      DEXSeqHTML( dxr, FDR=0.1, color=c("#FF000080", "#0000FF80") )

      每一个都是能点开的

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      (未测试)单细胞数据scanpy 学习笔记
      2020年2月23日

      下一篇文章

      (未测试)探索可变剪切分析6种方法
      2020年2月24日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?