> For the complete documentation index, see [llms.txt](https://zpliu.gitbook.io/booknote/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/di-7-ge-jie-guo.md). # 第7个结果 * 各个剪切事件并不是显著性的差异 * ES数据核实一下 * 同源基因之间的poly A的差异 * 组成型intron与剪切intron之间长度差异 * 同源基因在各个棉种的比较中AS的差异 * 同源基因间表达的差异与保守度 * 各个基因组中TE的插入 ## 处理A2的原始数据 参考 > AS数据的处理，由于D5基因组相比A2基因组，注释的转录本更多，这里看一下只考虑PacBio测序得到的与参考基因组进行比较的结果 > for i in A2 D5 TM-1; do for j in IntronR AltA AltD AltP ExonS Other; do printf ${j}"\t"; grep ${j} ../${i}/end\_third|grep PB|wc -l; done; done > > 和对应的基因数目 > > for i in A2 D5 TM-1; do for j in IntronR AltA AltD AltP ExonS Other; do printf ${j}"\t"; grep ${j} ../${i}/end\_third|grep PB|cut -f2|sort|uniq|wc -l; done; done | 基因组 | IR | AltA | AltD | AltP | ES | Other | | ---- | ----- | ---- | ---- | ---- | ---- | ----- | | A2 | 28078 | 5199 | 4135 | 1823 | 2276 | 1435 | | D5 | 22226 | 3257 | 2434 | 1131 | 1719 | 1039 | | TM-1 | 40354 | 7823 | 6854 | 4665 | 3071 | 2539 | | At | 19869 | 3868 | 3341 | 2301 | 1531 | 1197 | | Dt | 20436 | 3946 | 3505 | 2363 | 1540 | 1342 | ## 看看公司跑的数据，最后就使用这个数据了 > 鉴定可变剪切的原理: > > 1. 将PacBio测得的isoform比对到参考基因组，获得isform 注释信息 > 2. 将PacBio isoform信息与参考基因组组中原有的注释信息进行合并 > 3. 使用脚本，根据合并后的注释信息中的基因结构信息，提取对应剪切事件和发生的位置 > > 由于A2基因组中原有的注释信息里没有不同转录本的注释信息，并且PacBio测序只测了leaf这一个组织，所以在与A2中原有的注释信息合并后，也不能完全的代表所有的转录本的注释信息 > > 而D5和TM-1基因组中有不同转录本的注释信息，因此在根据基因结构去鉴定AS的时候，有有一些差异。 > > 为了减少这个带来的差异，我在各个棉种中选择那种只被PacBio isoform支持的AS > > 没有统计scaffold上的AS事件 * AS事件数目 | 基因组 | IR | AltA | AltD | AltP | ES | Other | | ---- | ----- | ---- | ---- | ---- | ---- | ----- | | A2 | 36574 | 5847 | 4472 | 4394 | 2352 | 4269 | | D5 | 30379 | 5350 | 3845 | 4185 | 2229 | 6302 | | TM-1 | 51453 | 9749 | 8252 | 7358 | 3723 | 7915 | | At | 25386 | 4878 | 4011 | 3621 | 1850 | 3882 | | Dt | 26067 | 4871 | 4241 | 3737 | 1873 | 4033 | * gene 数目 | 基因组 | IR | AltA | AltD | AltP | ES | Other | | --- | ----- | ---- | ---- | ---- | ---- | ----- | | A2 | 9139 | 3147 | 2614 | 2200 | 1499 | 1881 | | D5 | 7466 | 2900 | 2295 | 2214 | 1443 | 2285 | | TM1 | 13574 | 5414 | 4740 | 3780 | 2399 | 3308 | | At | 6657 | 2709 | 2345 | 1872 | 1210 | 1672 | | Dt | 6917 | 2705 | 2395 | 1908 | 1189 | 1636 | ### 提取FESTs序列 ```bash python ~/scripte/Alternative/module/extractFEST2.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 -AS ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt -g ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -IR 1 -ES 2 -AltA 3 -AltD 4 ## D5在分析的时候，取的是CDS坐标 ##拆分At与Dt ``` ### blast分析保守性 * 相似度80% ，覆盖度80%以上 e-value 1e-5 统计每种AS对应的比例 ```bash # evm.TU.Ga01G0013 0 Ghir_A01G000130 0 出现这种情况是因为，这个基因不存在IR事件，但是存在AltA获取AltD等其他事件 ## 统计同时存在AS的基因数，和AS数目 for i in `ls .`; do printf ${i}"\t"; awk '$2!=0&&$4!=0{a+=1;b+=$2+$3}END{print a"\t"b}' ${i}; done ## 统计保守的基因数，与AS数 ``` #### 保守IR事件 > 存在保守事件基因的比例： > > 再计算每个基因，AS事件的保守程度： > > 保守基因的比例: > > 1. 2425/7142 33.95% > 2. 2274/6594 34.48% > 3. 31.91% > 4. 62.59% > > 保守事件的比例: > > 1. 4789/36144 \*2 26.50% > 2. 26.95% > 3. 24.14% > 4. 50.14% | 比较 | 基因对数 | 事件数 | 保守基因对数 | 保守事件数没乘2 | | ---------- | ---- | ----- | ------ | -------- | | A2 vs At | 7142 | 36144 | 2425 | 4789 | | D5 vs Dt | 6594 | 32642 | 2274 | 4398 | | A2 vs D5 | 7218 | 35364 | 2303 | 4269 | | At vs Dt | 2911 | 13551 | 1822 | 3397 | #### 保守ES事件 > 保守基因的比例: > > 1. 31.10% > 2. 32.11% > 3. 23.28% > 4. 56.40% > > 保守事件的比例: > > 1. 33.85% > 2. 36.47% > 3. 26.42% > 4. 65.40% | 比较 | 基因对数 | 事件数 | 保守基因对数 | 保守事件数没乘2 | | ---------- | ---- | ---- | ------ | -------- | | A2 vs At | 598 | 1152 | 186 | 195 | | D5 vs Dt | 570 | 1146 | 183 | 209 | | A2 vs D5 | 524 | 984 | 122 | 130 | | At vs Dt | 250 | 474 | 141 | 155 | #### 保守AltA事件 > 保守基因的比例： > > 1. 55.22% > 2. 48.02% > 3. 41.05% > 4. 45.61% > > 保守事件的比例: > > 1. 48.47% > 2. 48.25% > 3. 41.08% > 4. 45.98% | 比较 | 基因对数 | 事件数 | 保守基因对数 | 保守事件数没乘2 | | ---------- | ---- | ---- | ------ | -------- | | A2 vs At | 891 | 2030 | 433 | 492 | | D5 vs Dt | 835 | 1919 | 401 | 463 | | A2 vs D5 | 765 | 1704 | 314 | 350 | | At vs Dt | 820 | 1792 | 374 | 412 | #### 保守AltD事件 > 保守基因的比例： > > 1. 49.92% > 2. 46.88% > 3. 40.26% > 4. 46.22% > > 保守事件的比例: > > 1. 52.32% > 2. 51.19% > 3. 44.55% > 4. 49.36% | 比较 | 基因对数 | 事件数 | 保守基因对数 | 保守事件数没乘2 | | ---------- | ---- | ---- | ------ | -------- | | A2 vs At | 659 | 1380 | 329 | 361 | | D5 vs Dt | 640 | 1344 | 300 | 344 | | A2 vs D5 | 544 | 1055 | 219 | 235 | | At vs Dt | 662 | 1414 | 306 | 349 | ### 比较不同的AS类型，保守程度是否存在差异 * 计算保守基因对，的AS事件保守程度 > - 不同剪切事件在多倍化的过程中，保守的程度不一样，其中ES最保守，而IR保守性最低 ![不同AS事件，保守性不一样](https://s1.ax1x.com/2020/05/14/YB0aD0.png) * 同一类剪切事件，在不同基因组的比较中的保守程度差异 > A2 vs At 叫保守，D5 vs Dt叫保守； > > A2 vs D5叫并行分化，At vs Dt叫趋同进化 ### 不同剪切事件的长度差异 * 组成型的内含子 * 保守的剪切事件对应的intronR * 不保守的剪切事件对应的intronR ```bash #所有的型内含子 extract_splice_sites.py ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf >1111 awk '{print $1,$2+2,$3,$4}' OFS="\t" 1111 >all_intron.bed awk '{print $1,$2+2,$3,$4}' OFS="\t" ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model_splce.txt >> all_intron.bed sort -k1,1 -k2,3n all_intron.bed |uniq >1111 mv 1111 all_intron.bed ~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ../all_intron.bed -b ./As_intronR.bed >222 # 发生剪切的内含子 grep ">" ../../../TM1/1|sed 's/>//g'|awk -F "-" '{print $2"\t"$3"\t"$4"\tAlternative"}' |uniq >As_intronR.bed #组成型的intron awk '$5=="."{print $0}$5!="."&&($7-$6)/($3-$2)<0.4{print $0}' 222|cut -f1,2,3|sort |uniq |awk '{print $0"\tConstitutive"}'>constitutive_intron.bed #存在剪切的内含子中保守的 cut -f3 ../../../blast/IR/A2_At_conserve_end|awk -F "-" '{print $2"\t"$3"\t"$4"\tconserve"}' >conserveAS.bed #存在剪切的内含子中不保守的 cat conserveAS.bed As_intronR.bed |sort |uniq -u >noconserveAS.bed ``` * 组成型的外显子 * 发生ES的外显子 * 保守的ES ```bash ##所有exon坐标 cat ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf |awk '$3~/exon/{print $1,$4,$5,$7}' OFS="\t" >all_exon.bed sort -k1,1 -k2,3 all_exon.bed |uniq >111 mv 111 all_exon.bed ## 发生剪切的exon的坐标 cat ../../../TM1/ES_length >./As.bed ##组成性exon cat all_exon.bed As.bed |cut -f1,2,3|sort |uniq -u |awk '{print $0"\tconstitutive"}' >constitutive.bed ## 保守的ES的坐标 cut -f3 ../../../blast/ES/A2_At_conserve_end|xargs -I {} grep {} As.bed |awk '{print $0"\tconserve"}' >conserve.bedd ## 不保守的ES坐标 cut -f1-4 conserve.bed|cat - As.bed |sort |uniq -u |awk '{print $0"\tnocnoserve"}' >noconserve.bed ``` ### 对多倍化过程中的剪切事件进行分类 * parallel conserve * A2与D5中保守，在At与Dt中同样保守 * A2与D5中不保守，但是在At与Dt中保守；A2 At或者D5 Dt中存在保守 * bias conserve * A2与D5中保守，但是在At与Dt中不保守；A2 At或者D5 Dt中存在保守 * A2与D5中不保守，在At与Dt中同样不保守；但是在A2 At或者D5 Dt中是保守的 ```bash ## 目录 ParallelEvolution ``` #### 四倍体内同源基因更多的发生分化 分析A2与D5中保守的剪切事件，在At、Dt中的保守情况 1. A2 D5中保守，在At与Dt中同样保守 2. A2 D5中保守，在At与Dt中不保守 A2和D5基因组分别是一个独立的个体，存在一些保守的AS，维持正常的生命活动； 剪切机制的存在往往是为了让植物能够快速响应外界环境的变化，在正常的环境中这种机制的存在往往是能量的浪费，因此在四倍体中，发现大多数原来在二倍体内保守的剪切事件，在四倍体内发生了分化，只在At中存在，或者只在Dt中存在，这两种比例各占一半 > A2和D5作为一个独立的个体，因此具有一些保守的剪切事件；而剪切事件往往是一种能量的浪费；在四倍体内只需要某一个基因组存在对应的AS即可 | 类型 | 二倍体内保守 | 四倍体内保守 | 四倍体内不保守 | | ------ | ------ | ------ | ------- | | IR基因 | 2303 | 536 | 1154 | | IR事件数 | 4269 | 838 | 1697 | | ES基因 | 122 | 21 | 46 | | ES事件 | 130 | 22 | 47 | | AltA基因 | 314 | 48 | 121 | | AltA事件 | 350 | 51 | 127 | | AltD基因 | 219 | 37 | 77 | | AltD事件 | 235 | 40 | 80 | > 绘制图形 > > 在画一个维恩图，A中保守的，D中保守的，At、Dt中同样保守的比例 ```bash ##A基因组中保守的事件数 cat 1 2 3 4 |grep "Ghir_A"|grep "evm"|awk '{print $1"^"$2"^"$3"^"$4}'>Agenome_conserve.txt cat 1 2 3 4 |grep "Ghir_D"|grep "Gor"|awk '{print $1"^"$2"^"$3"^"$4}' >Dgenome_conserve.txt ``` ![四倍体内存在分化](https://s1.ax1x.com/2020/05/17/YgcPED.png) ![事件数目统计](https://s1.ax1x.com/2020/05/18/Yh5x0g.png) ### 看剪切事件在两个亚基因组间是否存在偏向性 > 只在A基因组中保守的剪切事件，占A2和At中保守事件的比例 > > 只在D基因组中保守的剪切事件 | 事件类型 | 只在A保守的占A中所有保守比例 | 只在Dt中保守的占D中所有保守的比例 | | ---- | --------------- | ------------------ | | IR | 3623/4461 | 3189/4027 | | ES | 164/176 | 154/186 | | AltA | 394/445 | 358/409 | | AltD | 287/327 | 270/310 | ![剪切事件是否存在偏向性](https://s1.ax1x.com/2020/05/18/Yh57tA.png) ### 分析每个基因组中特异剪切事件对应的GO * 在A基因组中特异的AS对应的基因富集的GO * 在D基因组中特异的AS对应的基因富集的GO * A、D中保守的AS对应的基因富集的GO ```bash #在A基因组特异的 cut -f3 1|awk '{print "^^"$0"^"}'|cat - Agenome_conserve.txt | awk -F "^" '{print $3}'|sort|uniq -u|cut -f1 -d "-"|sort|uniq |xargs -I {} grep {} ~/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_Noscagenes_GO_V3_At.annot >A.GO #在D基因组中特异的 cut -f4 1 |awk '{print "^^^"$0}'|cat - Dgenome_conserve.txt |awk -F "^" '{print $4}'|sort|uniq -u|cut -f1 -d "-"|sort|uniq |xargs -I {} grep {} ~/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_Noscagenes_GO_V3_Dt.annot >D.GO ``` ### 分析可变剪切对基因表达式水平的影响 * At与Dt中存在保守事件的基因 * At与Dt中仅仅只有一个亚基因组中存在保守事件的基因 > 存在保守事件的同源基因对间的差异表达情况 > > 不存在保守事件的同源基因对的差异表达情况 > > Divergence of duplicated genes by repeated partitioning of splice forms and subcellular localization ```bash ## 获取同源基因对 cat 1 2 |awk '{print $3"\n"$4}' |awk -F "-" '{print $1}'|sort |uniq cat 3 4 |awk '{print $3"\n"$4}' |sed '/*/d'|awk -F "-" '{print $1}'|sort |uniq ## 获得同源基因对间的表达量 python homolog_expresion.py -homolog ~/work/Alternative/result/homologo/homologGene/At_vs_Dt_collinerity.txt -gene IR/At_Dt_bias_gene_ID -FPKM ~/work/RNA-seq/hisat2_out/leaf/leaf_BAM/Gh/gene_fpkm.txt -o IR/bias_fpkm ``` ## 剪切事件在多倍化过程中的变化 * 在二倍体中存在的事件，并且在四倍体中At、Dt保守；多倍化之前就存在的事件 * 在二倍体中不存在的事件，但是在四倍体中保守；多倍化后新产生的事件 ```bash # 上个分类中1、2 文件 # 新产生的做个排除法 ``` ## 参考 方差分析 卡方检验