可变剪切的鉴定

可变剪切的鉴定

为了方便在不同的基因组之间进行比较，AS的鉴定还是选择那些比较PacBio转录本之间的AS；从11_AS/end_splice.txt`文件中提取两个转录本都是PB注释的AS。将冗余的内含子保留信息合并一下

公司在鉴定AS的时候，把PB和参考基因组的转录本进行了合并；用参考基因组中原有的注释代替了PacBio的注释；而A2中只有一个转录本的注释。所以当我去取两个都是PB的转录本的时候，A2的就会比较多。

重新鉴定AS，使用gmap将PacBio比对到参考基因组输出格式设置为cDNA格式的gff文件

在使用AS的脚本去鉴定AS

  • test.gff 文件为Gmap的输出结果

  • PacBio.gff为PacBio转录本的注释信息

  • as 统计剪切位点处的编码

  • ats 输出AS鉴定的统计信息

  • os 输出剪切处的fasta序列

  • op 输出SVG图片信息

##先使用GMap进行比对
##再比对PacBio转录本之间的信息，使用脚本
/usr/bin/python ~/software/pipeline-for-isoseq/other/alternative_splice.py -i ./test.gff -g ./PacBio.gtf  -f ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -as -ats T -os  -t exon   -op -o /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/11_AS/alter_splice
##使用脚本对根据剪切编码进行分类
python ~/scripte/Alternative/classSplice.py
##根据去除冗余后的isoform提取其中的AS
python adjustAS.py -isoformId /public/home/zpliu/work/Alternative/result/homologo/FEST5/isoform/cdHit/A2_merge_isoform  -AS alter_splice2/111  -o test
###去除冗余
python ~/work/Alternative/result/homologo/FEST4/conserveAS/allAS/arrangeAS.py  -h
##统计数目

使用AS脚本鉴定AS

只鉴定PacBio转录本之间的IR，首先需要将../07_annotation/TM-1_merge_C.gtf文件中只提取PacBio的注释信息

##提取有基因注释的PacBio转录本信息，里面包含了lncRNA
grep "PB" ../07_annotation/TM-1_merge_C.gtf >PacBio.gtf
## 运行鉴定可变剪切的脚本

统计各个棉种中的数量信息，在没去冗余之前

基因组

IR

ES

AltA

AltD

AltP

Other

TM1

36172

2503

8033

6422

5876

4041

D5

22540

1323

4669

3617

3648

2465

A2

32179

2386

7186

5646

5553

3995

去完冗余转录本后，AS的统计

基因组

IR

ES

AltA

AltD

AltP

Other

TM1

33225

2208

6144

4848

3796

2549

D5

20680

1171

3497

2628

2343

1513

A2

29488

2120

5411

4256

3570

2448

使用rmats-turbo鉴定AS

安装软件

##使用conda环境
conda create  -n rmats-turbo python=3.6
conda activate rmats-turbo
##安装这些依赖项
pip install Cython==0.29.14 --user
pip install numpy==1.16.5 --user
conda install BLAS LAPACK -y
conda install GSL=2.5
#conda install GCC=5.4.0
conda install CMake=3.15.4
conda install gfortran

一行命令搞定

错误bash: gsl-config: command not found

主要是由于./build_rmats这个脚本运行过程中会source ~/.bashrc进入base环境就找不到脚本

解决办法把conda_envs/rmats/bin添加到环境变量里

装R包nloptr遇到一堆问题，最后还是用conda解决的

conda install -c r r-nloptr

##先加载GSL
module load GSL/2.4
##使用conda安装依赖
./build_rmats   --conda
## 安装R包
Rscript install_r_deps.R
##最后再编译一次
./build_rmats

使用

--statoff分析单个样本中的AS，不会进行多样吧间的比较

  • --b1指定比对·后的BAM文件所在位置 与--s1效果一样，只不过--s1需要调用STAR进行一次比对

  • --od 输出文件夹 outpu Dir

  • --tmp 临时文件所在目录

  • -t指定RNA-seq为单端或者双端测序

  • --libType 建库方式

  • --cstat对两个样本进行比较时的阈值

  • --nthread程序运行的线程数

##激活conda环境
conda activate conda_envs/rmats/
##运行脚本 
python ~/github/rmats-turbo/rmats.py  --b1 b1.txt  --gtf PacBio.gtf --readLength 150 -t paired --nthread 20 --od ./out --tmp ./tmp/ --statoff

输出文件

  • [AS_Event].MAT.JC.txt 跨越junction的read

  • [AS_Event].MAT.JCEC.txt 包括跨域junction的read和没有cross exon boundary的read

  • fromGTF.[AS_Event].txt 从GTF和RNA数据中鉴定到的AS

  • fromGTF.novelJunction.[AS_Event].txt再使用RNA的数据后才鉴定到的AS

每种AS类型中又有不同的列，最主要的差别就是指示事件的坐标，以及上下游外显子的坐标

具体的说明参考: https://github.com/Xinglab/rmats-turbo#output；

需要注意下面两个指标，在标准化的时候是以这两个为准，而不是真实事件的长度

https://groups.google.com/forum/#!topic/rmats-user-group/d7rzUBKXF1U

  • SkipFormLen这个只read支持跨域的实际长度，会和事件真实的坐标有差距

  • lncFormLen read cross整个坐标的长度

统计可变剪切数目

  • SE 外显子跳跃

  • 内含子保留

  • 可变的5‘

  • 可变的3’

  • 多个外显子互斥

###统计事件数目
for i in RI SE A5SS A3SS MXE; do 
awk -F "\t" '$4~/chrGhir/{print $0}' ${i}.MATS.JCEC.txt |wc -l; 
awk -F "\t" '$4~/chrGhir_A/{print $0}' ${i}.MATS.JCEC.txt |wc -l; 
awk -F "\t" '$4~/chrGhir_D/{print $0}' ${i}.MATS.JCEC.txt |wc -l;
done
##统计基因数目
for i in RI SE A5SS A3SS MXE; do 
awk -F "\t" '$4~/chrGhir/{print $0}' ${i}.MATS.JCEC.txt |cut -f2|sort|uniq |wc -l; 
awk -F "\t" '$4~/chrGhir_A/{print $0}' ${i}.MATS.JCEC.txt |cut -f2|sort|uniq  |wc -l; 
awk -F "\t" '$4~/chrGhir_D/{print $0}' ${i}.MATS.JCEC.txt |cut -f2|sort|uniq |wc -l;
done

这里的统计不包括scaffold

基因组

RI

SE

A5SS

A3SS

MXE

TM1

9735

4605

1839

2996

57

At

4793

2353

909

1488

27

Dt

4942

2252

930

1508

30

A2

7874

3707

1427

2391

56

D5

6262

2809

1172

1845

48

统计对应的基因数目

基因组

RI

SE

A5SS

A3SS

MXE

TM1

5562

3524

1506

2290

53

At

2736

1806

739

1141

25

Dt

2826

1718

767

1149

28

A2

4194

2638

1114

1767

50

D5

3559

2128

943

1376

41

使用阀值对AS进行一遍筛选

##外显子跳跃, 阀值小于0.9
awk -F "\t" '{split($(NF-2),a,",");if((a[1]+a[2])/2<=0.9){print $0}}' SE.MATS.JCEC.txt >filter/SE.JCEX.txt
awk -F "\t" '{split($(NF-2),a,",");if((a[1]+a[2])/2>=0.1){print $0}}' RI.MATS.JCEC.txt >filter/RI.JCEX.txt
awk -F "\t" '{split($(NF-2),a,",");if((a[1]+a[2])/2>=0.1){print $0}}' A5SS.MATS.JCEC.txt >filter/A5SS.JCEX.txt
awk -F "\t" '{split($(NF-2),a,",");if((a[1]+a[2])/2>=0.1){print $0}}' A3SS.MATS.JCEC.txt >filter/A3SS.JCEX.txt
awk -F "\t" '{split($(NF-2),a,",");if((a[1]+a[2])/2<=0.9){print $0}}' MXE.MATS.JCEC.txt >filter/MXE.JCEX.txt

high confidence AS event，筛选指标RI >=0.1 ES<=0.9

这里的统计不包括scaffold

基因组

RI

SE

A5SS

A3SS

MXE

TM1

6953

2135

1448

2489

54

At

3452

1094

732

1225

26

Dt

3501

1041

716

1264

28

A2

4891

1444

1149

1986

54

D5

4215

1076

916

1565

47

统计对应的基因数目

基因组

RI

SE

A5SS

A3SS

MXE

TM1

4175

1702

1198

1953

50

At

2068

875

607

964

24

Dt

2107

827

591

989

28

A2

2814

1086

909

1500

48

D5

2504

850

752

1191

40

参考

  1. STAR

  2. lr2rmats

  3. rmats-turbo

  4. rmats-turbo输出文件

Previousrmats2sashimiplotNext使用单个样本的数据进行AS分析

Last updated