明码标价之转录组测序数据的可变剪切

Original 生信技能树生信技能树 2022-08-10

收录于合集

#明码标价 24 个

#可变剪切 14 个

前面我们的明码标价之普通转录组上游分析，受到了各大热心粉丝的吐槽，觉得太简单了我们居然还好意思收费。

额，其实呢，这些粉丝应该是“饱汉不知饿汉饥”，其实数据分析这个技能啊，难者不会，会者不难。而且普通转录组上游分析仅仅是一个引子，又不是说我们只能提供这些粗浅的分析，转录组相关的绝大部分高级分析我们都是有现成的代码，而且都分享过一些教程：

用LeafCutter探索转录组数据的可变剪切，需要star软件的
使用SGSeq探索可变剪切，一个读取bam文件的R包

最近看到了一个比较新，而且发表的杂志也还不错的工具：SUPPA2

首先认识 SUPPA 软件，它有5个子命令，分别如下:

generateEvents: Generates events from an annotation.
psiPerEvent : Quantifies event inclusion levels (PSIs) from multiple samples.
psiPerIsoform : Quantifies isoform inclusion levels (PSIs) from multiple samples.
diffSplice : Calculate differential splicing across multiple conditions with replicates.
clusterEvents : Cluster events according to PSI values across conditions.

可以根据已知的 gtf文件信息针对多种可变剪切形式进行差异分析，它还提供了一些独立的py脚本文件，如下所示:

准备工作

包括conda安装，以及SUPPA2环境安装，参考转录组fasta序列获取和构建索引文件，参考gtf文件处理。

统一使用gencode最新版：https://www.gencodegenes.org/human/

mkdir -p SUPPA2/ref
cd SUPPA2/ref

nohup wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_37/gencode.v37.annotation.gtf.gz &
nohup wget  ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_37/gencode.v37.transcripts.fa.gz &
# 1.4G 2月  17 16:50 gencode.v37.annotation.gtf

使用conda创建rna小环境：

代码如下：

conda create -n rna
conda activate rna  
conda install -y trim-galore
# python-3.7.9         | 57.3 MB,  openjdk-10.0.2       | 189.2 MB
conda install -y star hisat2 bowtie2 # 没有其它依赖
conda install -y subread  bedtools deeptools
conda install -y salmon=1.4.0
conda install  -y -c bioconda suppa

下载参考转录组fasta序列文件然后使用salmon构建index文件

conda activate rna  
salmon index -t  gencode.v37.transcripts.fa  -i   gencode.v37.transcripts.salmon.index

得到的salmon的index是一个文件夹，内容如下:

$ ls -lh gencode.v37.transcripts.salmon.index/|cut -d" " -f 5-  

913K 2月  17 17:19 complete_ref_lens.bin
 95M 2月  17 17:23 ctable.bin
4.3M 2月  17 17:23 ctg_offsets.bin
175K 2月  17 17:19 duplicate_clusters.tsv
1.1K 2月  17 17:24 info.json
 90M 2月  17 17:24 mphf.bin
479M 2月  17 17:24 pos.bin
 496 2月  17 17:24 pre_indexing.log
 23M 2月  17 17:23 rank.bin
1.8M 2月  17 17:23 refAccumLengths.bin
 16K 2月  17 17:24 ref_indexing.log
913K 2月  17 17:23 reflengths.bin
 92M 2月  17 17:23 refseq.bin
 45M 2月  17 17:23 seq.bin
 126 2月  17 17:24 versionInfo.json

下载基因注释gtf文件然后使用suppa构建可变剪切事件综合文件

统一使用gencode最新版：https://www.gencodegenes.org/human/

对gtf文件的全部转录本走 suppa.py generateEvents 命令：

# 使用 generateEvents 命令根据基因组的gtf注释文件生成所有的可变剪切事件，格式保存为ioe格式
# 针对  SE SS MX RI FL   可变剪切批量分析的
suppa.py generateEvents \
-i   gencode.v37.annotation.gtf \
-o gencode.v37.events -e SE SS MX RI FL -f ioe

得到的文件如下:

$ wc -l *ioe
   20838 gencode.v37.events_A3_strict.ioe
   18921 gencode.v37.events_A5_strict.ioe
  112101 gencode.v37.events_AF_strict.ioe
   37317 gencode.v37.events_AL_strict.ioe
    8629 gencode.v37.events_MX_strict.ioe
    8474 gencode.v37.events_RI_strict.ioe
   50562 gencode.v37.events_SE_strict.ioe

可以看到是，基于 gencode.v37.annotation.gtf 文件的7种全部的可能的可变剪切事件，都被列出来了。

# 合并所有的ioe文件
awk 'FNR==1 && NR!=1 { while (/^<header>/) getline; } 1 {print}' *.ioe > gencode.v37.all.events.ioe
wc -l gencode.v37.all.events.ioe
# 256842 gencode.v37.all.events.ioe

上游分析拿到tpm矩阵

测序数据获取

得到文件如下:

ls -lh rawdata/*gz|cut -d" " -f 5-
1.7G 2月   8 14:33 rawdata/SRR1513329_1.fastq.gz
1.7G 2月   8 14:33 rawdata/SRR1513329_2.fastq.gz
1.7G 2月   8 14:33 rawdata/SRR1513330_1.fastq.gz
1.7G 2月   8 14:33 rawdata/SRR1513330_2.fastq.gz
1.8G 2月   8 14:33 rawdata/SRR1513331_1.fastq.gz
1.8G 2月   8 14:33 rawdata/SRR1513331_2.fastq.gz
1.8G 2月   8 14:33 rawdata/SRR1513332_1.fastq.gz
1.8G 2月   8 14:33 rawdata/SRR1513332_2.fastq.gz
1.9G 2月   8 14:33 rawdata/SRR1513333_1.fastq.gz
2.0G 2月   8 14:33 rawdata/SRR1513333_2.fastq.gz
1.7G 2月   8 14:33 rawdata/SRR1513334_1.fastq.gz
1.8G 2月   8 14:33 rawdata/SRR1513334_2.fastq.gz

然后质控，走fastqc软件，查看报告

fastqc -t 2 *.fastq.gz -o  ../qc

通常情况下大家需要针对clean后的fq文件进行后续的salmon流程定量获得各个样品的TPM值。

走trim流程

因为是双端测序，所以:

cd  rawData
 ls *gz|cut -d"_" -f 1|sort -u |while read id;do
nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ../cleanData   ${id}*.gz &
done

得到的clean的fq文件如下:

 $ ls -lh *gz|cut -d" " -f 5-
1.6G 2月  18 10:00 SRR1513329_1_val_1.fq.gz
1.6G 2月  18 10:00 SRR1513329_2_val_2.fq.gz
1.6G 2月  18 10:00 SRR1513330_1_val_1.fq.gz
1.6G 2月  18 10:00 SRR1513330_2_val_2.fq.gz
1.7G 2月  18 10:02 SRR1513331_1_val_1.fq.gz
1.8G 2月  18 10:02 SRR1513331_2_val_2.fq.gz
1.7G 2月  18 10:02 SRR1513332_1_val_1.fq.gz
1.7G 2月  18 10:02 SRR1513332_2_val_2.fq.gz
1.9G 2月  18 11:15 SRR1513333_1_val_1.fq.gz
1.9G 2月  18 11:15 SRR1513333_2_val_2.fq.gz
1.7G 2月  18 18:04 SRR1513334_1_val_1.fq.gz
1.7G 2月  18 18:04 SRR1513334_2_val_2.fq.gz

项目介绍分组情况如下:

negative control siRNA ：SRR1513329,SRR1513330,SRR1513331   
TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334

salmon流程定量获得各个样品的TPM值

## 定量获得TPM值  
mkdir salmon_output
cd cleanData  
index=../ref/gencode.v37.transcripts.salmon.index/
 ls *gz|cut -d"_" -f 1|sort -u |while read id;do
 nohup salmon quant -i $index  -l ISF --gcBias \
 -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -p 2 \
 -o  ../salmon_output/${id}_output 1>${id}_salmon.log 2>&1  &
 done
 ## quant.sf文件很重要，要用于后续的分析
##ENST和ENSG的前三个字母（ENS），意思是“ENSENMBLE”。T是指转录本。G是指基因。当然，也会看到ENSP，P自然是指蛋白质了。

几分钟就可以走完全部的salmon流程，每一个样品虽然说会输出一个文件夹，但是最重要的都是各自的 quant.sf文件，行数都是统一的:

$ wc -l   ../salmon_output/*/quant.sf
   233653 ../salmon_output/SRR1513329_output/quant.sf
   233653 ../salmon_output/SRR1513330_output/quant.sf
   233653 ../salmon_output/SRR1513331_output/quant.sf
   233653 ../salmon_output/SRR1513332_output/quant.sf
   233653 ../salmon_output/SRR1513333_output/quant.sf
   233653 ../salmon_output/SRR1513334_output/quant.sf

各自的 quant.sf文件内容如下:

==> ../salmon_output/SRR1513334_output/quant.sf <==
Name Length EffectiveLength TPM NumReads
ENST00000456328.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000362751.1|DDX11L1-202|DDX11L1|1657|processed_transcript| 1657 1587.157 0.000000 0.000
ENST00000450305.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000002844.2|DDX11L1-201|DDX11L1|632|transcribed_unprocessed_pseudogene| 632 452.000 0.000000 0.000
ENST00000488147.1|ENSG00000227232.5|OTTHUMG00000000958.1|OTTHUMT00000002839.1|WASH7P-201|WASH7P|1351|unprocessed_pseudogene| 1351 1273.856 3.098744 27.492

tpm表达量矩阵的获取

这个时候，可以直接使用作者的 multipleFieldSelection.py 脚本，对我们的全部样品的 quant.sf文件进行信息提取即可。其实自己写一个脚本也是很简单的：

-k indicates the row used as index
-f TPM值在第四列，在这里要提取每个样本的TPM值

## 提取所有样本的TPM值并合并为一个文件
 cd ../salmon_output/
multipleFieldSelection.py \
-i  SRR*/quant.sf -k 1 -f 4 \
-o iso_tpm.txt

得到 iso_tpm.txt的内容，还需要格式化一下，因为它跟gtf文件里面的基因信息不一致。

# gencode.v37.transcripts.fa 内容如下
>ENST00000450305.2|ENSG00000223972.5|OTTHUMG00000000961.2|OTTHUMT00000002844.2|DDX11L1-201|DDX11L1|632|transcribed_unprocessed_pseudogene|
# gencode.v37.annotation.gtf 内容如下
chr1    HAVANA  gene    11869   14409   .       +       .       gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; hgnc_id "HGNC:37102"; havana_gene "OTTHUMG00000000961.2";

各自的 quant.sf文件里面的ID是gencode.v37.transcripts.fa里面的转录本ID，有点过于复杂了，所以需要进行简化，需要与 gencode.v37.annotation.gtf 里面的转录本ID保持一致。

 perl -alne '{/(\|.*\|)\t/; ;s/$1//g;s/\|//g;print}' iso_tpm.txt > iso_tpm_formatted.txt

输出文件是 iso_tpm_formatted.txt ，是 233653 行的普通文本，可以看到简洁很多：

SRR1513329_output SRR1513330_output SRR1513331_output SRR1513332_output     SRR1513333_output SRR1513334_output
ENST00000456328.2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000450305.2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

拿到了这样的以转录本ID为行的tpm矩阵，再加上 gencode.v37.annotation.gtf 的可变剪切形式文件，就可以进行下游分析获取统计学显著的可变剪切事件。

下游分析获取统计学显著的可变剪切事件

项目介绍分组情况如下:

negative control siRNA ：SRR1513329,SRR1513330,SRR1513331   
TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334

对tpm表达量矩阵，根据可变剪切事件计算psi值矩阵：

ioe_merge_file=../ref/gencode.v37.all.events.ioe
ls $ioe_merge_file
suppa.py psiPerEvent \
-i $ioe_merge_file -e   iso_tpm_formatted.txt   -o project_events 1>psiPerEvent_log.txt 2>&1 &
# 会输出  project_events.psi 文件

同时可以检查看看 psiPerEvent_log.txt 文件。

project_events.psi 文件，内容如下所示：

SRR1513329_output SRR1513330_output SRR1513331_output SRR1513332_output     SRR1513333_output SRR1513334_output
ENSG00000000003.15;A5:chrX:100635746-100636191:100635746-100636608:- 0.20449107275352194   0.20506877462332357 0.18431225084229827 0.32174184236220055 0.27519806289263166   0.318119876931685
ENSG00000000003.15;A5:chrX:100635746-100636608:100635746-100636793:- 1.0 1.0 0.9009096685944438 1.0 0.9145937251391483 1.0

根据分组做差异可变剪切事件寻找

把前面的全部转录本可变剪切事件的psi文件以及 iso_tpm_formatted.txt ，根据分组信息，拆分为两个psi矩阵和两个tpm矩阵。

# negative control siRNA ：SRR1513329,SRR1513330,SRR1513331   
cut -f 1-4 project_events.psi > control.psi
cut -f 1-4 iso_tpm_formatted.txt > control.tpm
# TRA2A/B siRNA ：SRR1513332,SRR1513333,SRR1513334  
 cut -f 1,5-7 project_events.psi > treat.psi
  cut -f 1,5-7 iso_tpm_formatted.txt > treat.tpm
  # 很有意思的事情是，这个时候的表头似乎是并不是很重要

然后走 suppa.py diffSplice 命令，对两个分组的psi矩阵和tpm矩阵进行差异可变剪切事件寻找

ioe_merge_file=../ref/gencode.v37.all.events.ioe
ls $ioe_merge_file
suppa.py   diffSplice \
-m empirical -gc -i  $ioe_merge_file   \
--save_tpm_events \
-p treat.psi control.psi  \
-e treat.tpm control.tpm \
-o project_diffSplice # 前缀

差异分析生成下列3个文件：

ls -lh  project_diffSplice*    |cut -d" " -f 5-

7.9M 2月  18 09:16 project_diffSplice_avglogtpm.tab
 23M 2月  18 09:38 project_diffSplice.dpsi
 29M 2月  18 09:38 project_diffSplice.psivec

到此为止，统计学显著的可变剪切事件就获取完毕，后面就是各式各样的可视化啦！

当然了，一般来说，统计学显著的才有意义:

$ cat project_diffSplice.dpsi|perl -alne '{print if $F[2] <0.01}' |wc
    118     353   11552
$ cat project_diffSplice.dpsi|perl -alne '{print if $F[2] <0.05}' |wc
    750    2249   74204

查看差异比较大的可变剪切事件：

treat-control_dPSI treat-control_p-val
ENSG00000008517.18;AF:chr16:3065297:3065494-3065784:3065677:3065710-3065784:+ -1.0 0.0
ENSG00000009307.16;SE:chr1:114738072-114749821:114750207-114757925:- 0.6712142721 0.008991009
ENSG00000009307.16;SE:chr1:114739890-114749821:114750207-114757925:- 0.1488639077 0.008991009

这些差异比较大的可变剪切事件可以进行一些可视化，可以是igv里面的，举例如下:

网页工具绝对是完成不了这样的命令行数据分析哦

这个是基于Linux的ngs数据的上游处理，目前没有成熟的网页工具支持这样的分析。其实呢，如果你有时间请务必学习编程基础，自由自在的探索海量的公共数据，辅助你的科研，那么：

如果你没有时间从头开始学编程，也可以委托专业的团队付费拿到同样的数据分析，这样的20个样品以内的简单分组转录组测序数据的表达矩阵获取，首先需800元，然后如果要进行高级分析，比如上面讲解的可变剪切，再需要800元，仍然是可以拿到全部的数据和代码哦！

需要自己读文献筛选合适的数据集
提供半个小时左右的一对一讲解转录组数据处理背景知识。

如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。

李光耀：过早翘起尾巴与美国对抗是中国厄运的开始！

谁会想到，裁员会裁到总编辑头上

法官累积受贿929万：介绍案源、解冻账户、提取管理人报酬

现在有什么副业可以让人快速上岸? 可以试试这个行业，上岸其实不难!!

一场内斗，两败俱伤！

明码标价之转录组测序数据的可变剪切

准备工作

使用conda创建rna小环境：

下载参考转录组fasta序列文件然后使用salmon构建index文件

下载基因注释gtf文件然后使用suppa构建可变剪切事件综合文件

上游分析拿到tpm矩阵

测序数据获取

走trim流程

salmon流程定量获得各个样品的TPM值

tpm表达量矩阵的获取

下游分析获取统计学显著的可变剪切事件

根据分组做差异可变剪切事件寻找

网页工具绝对是完成不了这样的命令行数据分析哦

您可能也对以下帖子感兴趣

李光耀：过早翘起尾巴与美国对抗 是中国厄运的开始！

谁会想到，裁员会裁到总编辑头上

法官累积受贿929万：介绍案源、解冻账户、提取管理人报酬

现在有什么副业可以让人快速上岸? 可以试试这个行业，上岸其实不难!!

一场内斗，两败俱伤！

生成图片，分享到微信朋友圈

明码标价之转录组测序数据的可变剪切

准备工作

使用conda创建rna小环境：

下载参考转录组fasta序列文件然后使用salmon构建index文件

下载基因注释gtf文件然后使用suppa构建可变剪切事件综合文件

上游分析拿到tpm矩阵

测序数据获取

走trim流程

salmon流程定量获得各个样品的TPM值

tpm表达量矩阵的获取

下游分析获取统计学显著的可变剪切事件

根据分组做差异可变剪切事件寻找

网页工具绝对是完成不了这样的命令行数据分析哦

您可能也对以下帖子感兴趣

李光耀：过早翘起尾巴与美国对抗是中国厄运的开始！