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只提取包含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

结果部分

Previous02Next第三个结果

Last updated