> For the complete documentation index, see [llms.txt](https://zpliu.gitbook.io/booknote/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/di-wu-ci-fen-xi/ke-bian-jian-qie-de-jian-ding.md). # 可变剪切的鉴定 可变剪切的鉴定 为了方便在不同的基因组之间进行比较，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类型中又有不同的列，最主要的差别就是指示事件的坐标，以及上下游外显子的坐标 > 具体的说明参考: > > 需要注意下面两个指标，在标准化的时候是以这两个为准，而不是真实事件的长度 > > * `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/)