可变剪切的鉴定

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