# 可变剪切的鉴定

可变剪切的鉴定

为了方便在不同的基因组之间进行比较，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图片信息

```bash
##先使用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
```

```bash
###去除冗余
python ~/work/Alternative/result/homologo/FEST4/conserveAS/allAS/arrangeAS.py  -h
##统计数目
```

### 使用AS脚本鉴定AS

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

```bash
##提取有基因注释的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

### 安装软件

```bash
##使用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

```bash
##先加载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`程序运行的线程数

```bash
##激活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’
* 多个外显子互斥

```bash
###统计事件数目
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进行一遍筛选

```bash
##外显子跳跃, 阀值小于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](https://github.com/alexdobin/STAR)
2. [lr2rmats](https://github.com/Xinglab/lr2rmats)
3. [rmats-turbo](https://github.com/Xinglab/rmats-turbo/tree/v4.1.0)
4. [rmats-turbo输出文件](https://www.bioinfo-scrounger.com/archives/346/)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/di-wu-ci-fen-xi/ke-bian-jian-qie-de-jian-ding.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.