可变剪切统计

统计AS与计算PSI值

SUPPA

在ioe模式下生成的gtf文件用于在USC浏览器中加载；

conda activate suppa
##提取感兴趣的AS的坐标
suppa.py  generateEvents -i PacBio.gtf -o TM1 -e SE SS RI FL MX -f ioe

D5总共鉴定到16043个剪切事件

A2总共鉴定到23245个剪切事件，去除冗余后对应了7725个基因

TM1总共鉴定到27595个剪切事件，去除冗余 后对应了10042个基因

基因组	RI	SE	A5	A3	Other	Total
D5	9981	2605	4702	8232	3549	29069
A2	12191	1964	4688	6673	2986	28502
TM1	16472	3111	7297	10418	3763	41061
At	8072	1562	3626	5158	1963	20381
Dt	8400	1549	3671	5260	1800	20680

38644

对应的基因数目

基因组	RI	SE	A5	A3	Other	Total
D5	4915	2047	3283	5209	1109
A2	5763	1486	3139	4159	941
TM1	8160	2446	4984	6821	1232	23643
At
Dt

对应的AS转录本数目

for i in `ls .|grep -E D52_.*ioe`; do awk '$4~/PB/{print $4}'  ${i} |sed 's/\,/\n/g'|sort |uniq |grep PB|wc -l; done

基因组	RI	SE	A5	A3	Other
TM1	17193	5571	12232	17112	2493
A2	13163	3762	8390	12246	2218
D5	10147	5140	7847	13787	2346

分析两种转录本在不同剪切事件的比例

提取每种剪切事件对应的转录本ID号在对它的性质进行判断，判断它是否是Annotion、unAnnotion的

##提取对应事件的PacBio编号，并且打上Annotion标签
for i in RI SE AL AF A3 A5 MX; do awk '$4~/PB/{print $4}' TM2_${i}_strict.ioe |sed 's/,/\n/g'|grep PB|awk '{print "'${i}'""\t"$0}'|sort |uniq  >${i}_isoform.txt; python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py  ../collapse/PacBio_Annotion/all_isoform.txt ${i}_isoform.txt 11; mv 11 ${i}_isoform.txt; done
##统计数目
for i in `ls .|grep isoform`; do awk '$(NF)=="unAnnotion"{a+=1}$(NF)=="Annotion"{b+=1}END{print a"\t"b}'  ${i}; done

分析发生AS的转录本与没发生AS的转录本表达量的差异

##给转录本打上AS的标签
cat *isoform.txt |awk '{print $0"\tAS\t"$2}'|sort -k5,5|uniq -f4 >Alternative_isoform.txt
##分析表达水平的差异
python ../ORF/AddAnnotionTag.py  Alternative_isoform.txt  ../collapse/PacBio_Annotion/all_isoform_FPKM all_isoform_FPKM
##统计AS的转录本的平均表达水平
awk '$(NF)=="AS"{print $(NF-2)}' ../AS2/all_isoform_FPKM |awk '{a+=$0}END{print a/NR}'