01

进一步对可变剪切进行分类

由于只测了叶片一个组织，加上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
## 统计每种事件的数目

棉种	IR	ES	AltA	AltD
A2	35756	1617	4242	3323
D5	26185	1616	3595	2633
TM1	43581	2455	6333	5176
At	21540	1199	3150	2555
Dt	22041	1256	3183	2621

发生AS的isform数目

基因组	isform数
A2	17538
D5	11812
At	10783
Dt	11171

检测到的转录出的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

基因组	isform数
A2	68833
D5	53097
TM1	83087
At	40789
Dt	42167

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

PacBio鉴定到的isform剪切起始位点与基因边界进行比较

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数

对比对带参考基因组上的isform 进行校正

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与基因的坐标

• 使用校正后的isform与gene的关系

• 如果isform坐标在gene坐标内，则保留

• 如果超出，则只保留到那个有交集的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与基因坐标比较

根据校正情况的转录本的表达水平进行校正

首先计算不同转录本的表达量

对注释文件中的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

鉴定保守的isform

每条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数目上没有明显的差异。

不同比较	isfrom数目	isform数目
A2 vs D5	55764	43577
At vs Dt	35529	36050
A2 vs At	57144	35182
D5 vs Dt	47685	38060

校正后的结果

不同比较	isfrom数目	isform数目
A2 vs D5	55764	43577
At vs Dt	35529	36050
A2 vs At	57144	35182
D5 vs Dt	47685	38060

对isform数目存在差异的基因进行GO富集分析

鉴定不同基因组中同源基因间保守的isform

FSM（full splice match）

• 保守的isform

  exon数目要相同、长度允许几个碱基的误差；第一个exon和最后一个exon的长度不进行考虑

• 特异性的isform

##获取同源基因间保守的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降级

统计不同基因组间shared isform数与特有的FSM isform数

• FSM isform

  这些isform内部的exon剪切是一样的，但是两端的exon会存在不同的TTS和ployA；就可以分析不同的TTS和ployA对基因转录的影响

  不同比较	第一个总isform数	第二个总的isfrom数	第一个基因数	第二个基因数	第1个匹配的FSM数	第二个匹配FSM数	第一个对应的基因数	第二个对应的基因数
  At Dt	30498	30484	10492	10492	9296	9291	5549	5549
  A2 D5	50106	40833	13596	13596	16486	15476	8714	8714
  A2 At	47638	33283	12154	12154	15619	13302	7773	7773
  D5 Dt	40980	35224	12705	12705	14918	13594	8058	8058

• 独有的isform 、包括另外一个基因组中不存在的isfrom的情况

• 只存在单个isform的基因中保守的频率和多个isform保守的频率;基因转录出isform的数目与FSM的数目

每个基因对应的isform的剪切率

1. 首先根据找到发生AS事件的isform，这些isform是叶片组织中经过AS产生的isform

2. 在不同基因组之间，比较发生AS的isform的保守性，并且比较它们的剪切活性

   • 首先isform是经过AS产生的

   • 对应的AS事件是保守的，则认为isform也是保守的

   • 如果由AS产生的isform不保守，则认为这个isform是棉种中特异的AS产生的isform

结果部分

• At Dt 合并取uniq后和祖先进行比较

• At和Dt中序列保守中isform的差异，不保守的序列差异

• 可变剪切与基因家族分析

• https://blast.advbiocomp.com/licensing/personal.html)