进一步对可变剪切进行分类
由于只测了叶片一个组织，加上A2基因组注释信息的不完整；现在对于AS的鉴定只保留发生在PacBio isform上的剪切事件，它也是代表了叶片组织中剪切事件的信息
python ~/scripte/Alternative/module/FEST3/ifPacBioAS.py -p ~/work/Alternative/result/Gh_result/CO31_32_result/06_Alignment/all.collapsed.gtf -r ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf -o 1111 -AS end_splice.txt
## D5的时候对参考基因组的注释问文件做了修改
##统计发生在PacBio isform上的AS
awk '$6~/PB/&&$8==6{print $0}$7~/PB/&&$8==7{print $0}' 1111
## 统计每种事件的数目
发生AS的isform数目
检测到的转录出的isform数目，不包括预测的新基因
##07文件夹中的merge.gtf
awk '$3~/transcript/&&$2~/PacBio/{print $0}' merge.gtf |grep "orginal_gene_id \"P" -v |wc -l
## 统计每个基因所表达出的isform
python ~/scripte/Alternative/module/FEST3/staticgeneIsform.py -gtf ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -o gene_isformCount.txt
PacBio鉴定到的剪切位点与参考基因组注释的剪切位点进行比较
## 提取PacBio isform gtf文件
grep "orginal_gene_id \"P" -v ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf >PacBio_isform.gtf
## 提取剪切位点
extract_splice_sites.py PacBio_isform.gtf >PacBio_splice.txt
## 提取二代测序支持的剪切位点数
cat PacBio_splice.txt ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model_splce.txt|sort -k1,1 -k2,3n|uniq -c|awk '$1==2{print $0}' |wc -l
python ~/scripte/Alternative/module/FEST3/isformcompareGeneBoundary.py -PacGtf ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -refGff ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gff3 -o isfrom_refgene.txt
## 统计整个isform完全在gene内的数目
awk '$2>=$6&&$3<=$7{print $0}' isfrom_refgene.txt |wc -l
##转录起始位点超出 gene边界的数目
awk '$4=="+"&&$2<$6||$4=="-"&&$3>$7{print $0}' isfrom_refgene.txt
## 转录终止位点超出gene边界的数目
awk '$4=="+"&&$3>$7||$4=="-"&&$2<$6{print $0}' isfrom_refgene.txt |wc -l
## 计算每个isform的二维坐标
awk '$4=="+"{print $1"\t"1-($7-$2+1)/($7-$6+1)"\t"($3-$6+1)/($7-$6+1)}$4=="-"{print $1"\t"1-($3-$6+1)/($7-$6+1)"\t"($7-$2+1)/($7-$6+1)}' isfrom_refgene.txt
## 超出范围的gene boundary范围的到时候，用二代数据看看有没有read覆盖
## 超出范围的isform是否是比对到了非gene区域
~/software/bedtools2-2.29.0/bin/subtractBed -a chrosomoe.bed -b ref_geneBody.bed >nonGeneBody.bed
## 提取PacBio isform 转录本的bed坐标
awk -F "\t" '$2~/PacBio/&&$3~/transcript/{print $1"\t"$4"\t"$5"\t"$9}' PacBio_isform.gtf |awk -F ";" '{print $1"\t"$2}'|awk -F "\t" '{print $1"\t"$2"\t"$3"\t"$5}'|sed -e 's/\"//g' -e 's/ transcript_id //g'
##看PacBio与非gene区域的交集情况
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a PacBio_Transcript.bed -b nonGeneBody.bed >1111
## 转录本有部分比对到了另外一个gene区域，删除;
awk '$5!="."{print $0}' 1111 |awk '$2<$6&&$3>$7{print $0}' |cut -f4|sort |uniq >delete_isformt.txt
## isform超出基因范围，但是超出的部分只在非gene区域，没有跨域另外一个基因
awk '$5!="."{print $0}' 1111 |awk '$2>=$6&&$3>=$7||$2<=$6&&$3<=$7{print $0}'|cut -f4|sort |uniq >777
cat delete_isformt.txt 777|sort |uniq >888
cat 888 delete_isformt.txt |sort |uniq -u >999
## 注释文件中基因与基因之间重叠导致没检查出的isform也删除
awk '$2<0||$3>1{print $1}' 222 |sort |uniq |cat - 888 |sort|uniq -u >>delete_isformt.txt
##对比对到intergenetic区域的isform使用二代数据进行检查
#提取对应的intergenetic 区域
cat 999 |xargs -I {} grep -E "{}\s+" 1111 |awk '{print $5"\t"$6"\t"$7"\t"$4}' >test.bed
# 与二代数据取交集
~/software/bedtools2-2.29.0/bin/intersectBed -a test.bed -loj -b ~/work/RNA-seq/hisat2_out/leaf/leaf_BAM/SRR5886147sorted.bam_filter.bed >intersect.txt
## 统计每个isform比对上的read数
python ~/scripte/Alternative/module/FEST3/correctionPacBio.py -PacBio1 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -ref ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gff3 -o tmppp
根据校正的情况，对可变剪切进行校正
发生剪切的isform属于同一个基因，先对gene，进行校正
python ~/scripte/Alternative/module/FEST3/correctionAS.py -isform gene_isformCount.txt2 -ref ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gff3 -AS ~/work/Alternative/result/Gr_result/CO41_42_result/11_AS/end_splice.txt2 -o 5555
提取isfrom与基因的坐标
如果超出，则只保留到那个有交集的exon坐标或者保留那个完全超出的一个exon。
python ~/scripte/Alternative/module/FEST3/isformcompareGeneBoundary.py -PacGtf ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -refGff ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gff3 -isform gene_isformCount.txt2 -o 111111
##对isform进行坐标化
awk '$4=="+"{print $1"\t"1-($7-$2+1)/($7-$6+1)"\t"($3-$6+1)/($7-$6+1)}$4=="-"{print $1"\t"1-($3-$6+1)/($7-$6+1)"\t"($7-$2+1)/($7-$6+1)}' 111111
##绘制点图
根据校正情况的转录本的表达水平进行校正
首先计算不同转录本的表达量
对注释文件中的isform只提取包含PacBio isform的注释
## isform.gtf文件还需要进一步处理，将基因名换成参考基因组的
awk '$2=="PacBio"{print $0}' ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf |awk -F "\t" '{split($9,a,";");print $1,$2,$3,$4,$5,$6,$7,$8,a[3]"; "a[2]";"}' OFS="\t" |sed 's/orginal_//g'>isform.gtf
stringtie ~/work/RNA-seq/hisat2_out/leaf/leaf_BAM//SRR5886147sorted.bam_filter --rf -G isform.gtf -p 10 -b ./
针对isform与gene的对应关系提取isform的表达量数据
python ~/scripte/Alternative/module/FEST3/correctionGeneExpression.py -isform ../gene_isformCount.txt2 -fpkm t_data.ctab -o isform_fpkm.txt
每条PacBio isform中被FL read支持的read数
CO31_32_result/06_Alignment/../06_Alignment/all.collapsed.rep.fa
至少一条FL read与isform有关，则认为它表达了
#fasta文件
>PB.1.1|Ghir_A10:229-4314(+)|c53146/f2p0/926 c53146/f2p0/926 isoform=c53146;full_length_coverage=2;non_full_length_coverage=0;isoform_length=926
grep ">" all.collapsed.rep.fa |awk -F "|" '{print $1"\t"$3}'|awk -F ";" '{print $1"\t"$2"\t"$3}'|awk -F "\t" '{print $1"\t"substr($3,22)"\t"substr($4,26)}'|sed 's/>//g'
##得到结果：全长FLreads支持的数目 非全长NFL read支持数目
PB.1.1 2 0
PB.1.2 1 0
PB.2.1 1 5
亚基因组同源基因间isform数目统计
##统计同源基因间表达出的isform数目和编号
python ~/scripte/Alternative/module/FEST3/statichomologgeneisform.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt -isform1 ../A2/gene_isformCount.txt -isform2 ../TM1/gene_isformCount.txt -o 111
sort 111 |uniq >gene_isform_count.txt
### 统计所有同源基因所表达出的isform数目
cut -f2 gene_isform_count.txt |sed '/^0$/d'|awk -F "," '{a+=NF}END{print a}'
可以发现在At、Dt两个亚基因组合并到同一个细胞核后，与祖先基因组相比，每个基因产生的isform数目发生了都发生了减少，很大程度上是由于At、Dt亚组间功能冗余只需要某一个亚基因组产生对应的isform即可
两个亚基因组单独存在是isfrom的数目同样存在很大的差异，但是在At、Dt聚合到同一个细胞核后基因转录的isform数目上没有明显的差异。
校正后的结果
FSM（full splice match）
保守的isform
exon数目要相同、长度允许几个碱基的误差；第一个exon和最后一个exon的长度不进行考虑
##获取同源基因间保守的isform
python ~/scripte/Alternative/module/FEST3/conserveHomologIsform.py -homolog gene_isform_count.txt -isform1 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -isform2 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -error 5 -o 2222
转录起始位置相差的碱基数 ploy相差的碱基数
Ghir_A01G000100 PB.5351.1 Ghir_D01G000120 PB.21439.1 40 -15
Ghir_A01G000230 PB.5357.1 Ghir_D01G000240 PB.21446.2 -5 10
## 随机获取exon序列
python ~/scripte/Alternative/module/FEST3/validateconserveIsform.py -con homolog_FSM_isform.txt -PacBio1 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -PacBio2 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -fa1 ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -fa2 ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -o blast.fa -o2 random_validate_isform.txt
随机抽取1000个基因对应的FSM isform对进行blast比对，发现多少比例的保守
> >
At vs Dt : 987/1000
A2 vs At: 997/1000
D5 vs Dt 1000/1000
A2 vs D5 979/1000
## 根据得到的序列进行blast
blastn -query blast.fa -db ./blast -outfmt '6 qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore' -out blast.out
## 相似度95以上，覆盖度95以上；获取FSM isform对
awk '$8>=95&&$9>=95{print $0}' blast.out |awk '$1~/Ghir_A/&&$2~/Ghir_D/{print $1"\t"$2}'|sort|uniq >1
## 看存在随机抽取的isform中有多少是真正的FSM isform对
cat random_validate_isform.txt |xargs -I {} grep {} 1 |wc -l
同一个基因组内的FSM
计算由同一个基因产生的FSM数，和对应isform的FPKM值
##使用脚本得到同一个基因的FSM isform
##获取对应的isform编号
awk '$1==$3&&$2!=$4{print $0}' 2222
##获得两两FSM比较中TTS和ployA的绝对距离
awk '$1==$3&&$2!=$4{print $1"\t"$2"\t"$3"\t"$4"\t"sqrt($5*$5)"\t"sqrt($6*$6)}' 2222 |cut -f5,6|sort|uniq |head
分析存在表达的isform 中nonsense-mediated mRNA decay 的比例与不存在表达的isform的无意义介导的降解
#分析这些FSM中的表达情况与无意义介导mRNA降级
FSM isform
这些isform内部的exon剪切是一样的，但是两端的exon会存在不同的TTS和ployA；就可以分析不同的TTS和ployA对基因转录的影响
独有的isform 、包括另外一个基因组中不存在的isfrom的情况
只存在单个isform的基因中保守的频率和多个isform保守的频率;基因转录出isform的数目与FSM的数目
首先根据找到发生AS事件的isform，这些isform是叶片组织中经过AS产生的isform
在不同基因组之间，比较发生AS的isform的保守性，并且比较它们的剪切活性
对应的AS事件是保守的，则认为isform也是保守的
如果由AS产生的isform不保守，则认为这个isform是棉种中特异的AS产生的isform
结果部分
At和Dt中序列保守中isform的差异，不保守的序列差异
