> 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-liu-ci-fen-xi/quan-chang-zhuan-lu-ben-shu-ju-de-tong-ji.md).

# 全长转录本数据的统计

## 全长转录本数据统计

> <https://www.cnblogs.com/jinhh/p/8328818.html>

​ 对于每个棉种的原始下机数据进行统计，统计raw sequencing data中read数和consensus read数

```bash
 samtools view m54136_180615_020020.subreads.bam |cut -f1|awk -F "/" '{print $NF}'|awk -F "_" '{print $2-$1+1}'|awk '{a+=$1}END{print a}'
```

| 棉种    | raw sequence base | consensus reads count |
| ----- | ----------------- | --------------------- |
| Ga\_1 | 12083371005       | 320874                |
| Ga\_2 | 12402699575       | 292447                |
| Gr\_1 | 9756254646        | 238566                |
| Gr\_2 | 9939418624        | 248497                |
| Gh\_1 | 16310848287       | 319380                |
| Gh\_2 | 15711398953       | 366003                |
| Total | 76GB              |                       |

1.
2. **原始数据到最终得到的转录本数**

   > 这个数据在原始数据的统计表1中有
3. 饱和曲线分析，表明Iso-seq在基因的检测上已经达到了饱和，而对于新转录本的检测仍旧还需要提高测序深度
4. **PacBio转录本与reference转录本长度的比较**
5. **PacBio转录本的exon数目与reference exons数目比较**
6. **比对转录本进行分类**

   > 这里可以用实验去验证一些参考基因组中没有注释的剪切位点的信息；
   >
   > 同时说明PacBio可以检测到一些新的剪切位点的信息
   >
   > 转录本的长度分布，最长的转录本和最短的转录本对应的基因：

   | 基因组 | incomplete-splice matches | unAnnotion isoform | None gene region | total | prfect  match |
   | --- | ------------------------- | ------------------ | ---------------- | ----- | ------------- |
   | TM1 | 46127                     | 37265              | 4924             | 88316 | 29737         |
   | A2  | 28544                     | 38570              | 1880             | 68994 | 15035         |
   | D5  | 34829                     | 17141              | 3264             | 55234 | 22787         |

   将Annotion进行细分：

   * 与某一个转录本的剪切位点完全相同
   * 与参考转录本的剪接位点部分相同

| 基因组 | 最短Isoform                         | 最长Isoform                           |
| --- | --------------------------------- | ----------------------------------- |
| A2  | evm.TU.Ga02G1745 302bp            | evm.TU.Ga10G0940 12893bp unAnnotion |
| D5  | Gorai.007G076700 301bp unAnnotion | Gorai.002G214500 12207bp            |
| TM1 | Ghir\_D13G001640 311bp            | Ghir\_A07G019100 15431bp unAnnotion |

> 测得的转录本长度范围在301bp到1543bp，对应的基因注释为
>
> * AT3G01490 301bp 蛋白激酶超家族
> * AT3G02260 15431bp 生长素转运蛋白

### 基因所转录出的PacBio转录本数

1. **对基因的转录本类型进行分析**

   * 只转录出Annotion 转录本的基因

   * 只转录出unAnnotion转录本基因

   * 转录出Annotion与unAnnotion转录本基因

   > 通过PacBio测序，发现许多基因中仍旧可以鉴定到一些没有被注释的剪切信息，

   | 基因组 | 只转录出Annotion 的 | 只转录出unAnnotion | 转录出Annotion与unAnnotion | total |
   | --- | -------------- | -------------- | ---------------------- | ----- |
   | TM1 | 14953          | 7861           | 9153                   | 31967 |
   | A2  | 6857           | 6200           | 7156                   | 20213 |
   | D5  | 9470           | 2137           | 6630                   | 18237 |

   > 平均每个基因被测到FL-read数目，和转录本数

   | 基因组 | 平均转录出Isoform数 | 平均测到FL-read数 |
   | --- | ------------- | :----------: |
   | A2  | 2.6           |     20.2     |
   | D5  | 3.3           |     16.5     |
   | TM1 | 2.8           |     14.8     |

#### 比较不同基因组中在单外显子，和多外显子基因上PacBio转录本类型的数目。

基因是否是单外显子还是多外显子，取决于参考基因组中exon数目最多的那个转录本。

> 多外显子基因转录出更多类型的PacBio转录本，同时单外显子基因平均转录PacBio转录本的数目大于1，说明基因在转录过程中存在丰富的剪接形式。
>
> 与此同时两个二倍体物种相比于四倍体物种有着更加多的转录本

```bash
##统计每个基因的exon数目
awk  -F "\t" '$3~/exon/{print $9}' ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf|sed -e 's/;.*//g' -e 's/transcript_id \"//g' -e 's/\"//g'|sort |uniq -c |sed 's/\..*//g'|sort -k2,2 -k1,1nr |uniq -f1 1 |awk '{print $2"\t"$1}' >gene_exon_count.txt
##统计每种exon基因转录出PacBio转录本的数目
##画图的时候去除掉没有检测到PacBio转录本的基因数
```

## 分析PloyA位点和剪接位点附近的motif

### 剪接位点2bp的序列

> 对于鉴定到的所有外显子坐标进行motif的分析，第一个核苷酸最后一个核苷酸位点附近的2bp序列作为受体位点和供体位点。分析内含子剪接位点的motif
>
> * 供体5‘ 端GU   DNA代码`GT`
> * 受体3’端 AG  DNA代码 `AG`
> * 多少内含子是规范的GU-AG

```bash
##提取PacBio转录本exon的5'和3’处的2bp序列
awk '$2~/PacBio/{print $0}' ../../CO31_32_result/07_annotation/merge.gtf >PacBio.gtf
##提起内含子边界坐标序列
extract_splice_sites.py PacBio.gtf >111
awk '$4=="+"{OFS="\t";print $1,$2+1,$2+3,$1"_"$2"_"$3"*5","1",$4;print $1,$3-2,$3,$1"_"$2"_"$3"*3","1",$4}$4=="-"{OFS="\t";print $1,$3-2,$3,$1"_"$2"_"$3"*5","1",$4;print $1,$2+1,$2+3,$1"_"$2"_"$3"*3","1",$4}'  111 >intron_boundary.bed 
#奇数行为5‘ 偶数行为3’
fastaFromBed  -bed intron_boundary.bed  -fi genome.fa -name+ -s -fo  intron_boundary.fa
##供体端的motif 比例
grep "*5::" -A1  intron_boundary.fa |sed '/^[>-]/d'|awk '$1=="GT"{a+=1}END{print a/NR}'
##受体端的motif 比例
grep "*3::" -A1  intron_boundary.fa |sed '/^[>-]/d'|awk '$1=="AG"{a+=1}END{print a/NR}'
##intron 规范边界比例
sed '/^>/d' intron_boundary.fa |awk 'NR%2!=0{printf $1}NR%2==0{print $0}'|awk '$0=="GTAG"{a+=1}END{print a/NR}'
```

| 基因组 | 规范供体 | 规范受体  | 规范的内含子 |
| --- | ---- | ----- | ------ |
| TM1 | 0.93 | 0.95  | 0.92   |
| A2  | 0.93 | 0.945 | 0.92   |
| D5  | 0.92 | 0.934 | 0.909  |

> Iso-seq 精确的确定了ployA的位置和内含子的边界
>
> 内含子边界的motif序列信息，以及主要的内含子是由哪种motif组成的。
>
> * 单外显子转录本
> * **多外显子基因，**&#x5176;中多少基因会产生多个转录本结构
> * 分析多外显子转录本的长度分布
> * 分析参考基因组中鉴定到的多个TSS基因在PacBio有多少重叠。

* `R1`区域外显子上游的内含子35bp序列
* `R2`区域外显子上游和下游的32bp序列
* R3区域外显子下游的内含子区域40bp序列
* PAS（每条转录本的最后的位置）上游的35bp用于搜索ployA信号。

**分析基因是否存在多TTS和多ployA基因**

> 分析存在多TTS和多ployA的基因，以及ployA的motif分类。
>
> 每个转录本的截取CS（切割位点）上游50bp的序列，使用`SignalSleuth2`，扫描CS位点上游1-40bp的motif序列
>
> 分析物种中特异性的ployA位点。

多倍化过程中同源基因的ploy或者TSS数目发生变化

> ploy信号到时候就将提取上游50bp序列，扫描上游40bp的序列
>
> top 10的motif序列，规范的序列 AATAAA ATTAAA

```bash
##提取CS位点上游的50bp序列
python ~/github/zpliuCode/Isoseq3/02spliceSitemotif/polyadenylation.py  ~/work/Alternative/result/Gh_result/CO31_32_result/collapse/merge.gtf ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta TM1_polyadenylation_50b.txt 
##扫描CS位点上游1-40bp的motif
+ -k 6 对6聚体motif进行排名
+ -once motif只计数一次
+ -topn 输出前10个频率最高的motif
perl ~/software/SignalSleuth2/SignalSleuth2.pl -seqfile TM1_polyadenylation_50b.txt -from 10 -to 50 -once T -sort T -topn 10 -cnt T -pssm T -k 6
统计出现次数
awk '{print $1"\t"$NF}' *.cnt
```

| motif排名 | TM1  | A2   | D5   |
| ------- | ---- | ---- | ---- |
| ATTTTT  | 9246 | 7270 | 5613 |
| TATTTT  | 8785 | 6848 | 5348 |
| AATTTT  | 8606 | 6844 | 5362 |
| TTTTAT  | 8504 | 6683 | 5305 |
| TTATTT  | 8372 | 6516 | 5077 |
| AATAAA  | 8044 | 6512 | 4800 |
| AAATTT  | 8016 | 6511 | 4893 |
| TTTATT  | 8014 | 6279 | 4945 |
| TTTTTT  | 7663 | 6001 | 4735 |
| ATATTT  | 7643 | 6114 | 4488 |

> 在TM-1、A2、D5中分别检测到8044、6512、4800个规范的ployadentlation motif AATAA

分析基因包含多个`ployA`和多个`TSS`位点的信息

在TM1、A2、D5中比较基因转录出转录本的个数：

> 58%,70%,65%的基因存在多个转录本，这些存在多个转录本的基因中，存在多个ployA或者TTS位点的基因。
>
> * 内部的剪接位点完全相同，存在只由于ployA或者TSS造成转录本的差异的比例。

| 基因组 | 单个转录本的基因 | 多个转录本的基因 | 基因数   | 转录本数  |
| --- | -------- | -------- | ----- | ----- |
| TM1 | 13361    | 18606    | 31967 | 83392 |
| A2  | 6046     | 14167    | 20213 | 67113 |
| D5  | 6326     | 11911    | 18237 | 51964 |

```bash
##比较基因内转录本间的剪接位点和TSS、ployA
python  ~/github/zpliuCode/Isoseq3/05polyadenylation/transcript_different_polyadenylation.py  ~/work/Alternative/result/Gr_result/CO41_42_result/collapse/merge.gtf D5.txt
##单个ployA的基因数
awk '$7==1{print $0}' TM1.txt |cut -f1|sort |uniq |wc -l
##单个TSS的基因数
awk '$8==1{print $0}' TM1.txt |cut -f1|sort |uniq |wc -l
##剪接位点相同的转录本数
awk '$6=="splitsame"{print $2"\n"$3}' TM1.txt |sort |uniq |wc -l
```

**多转录本的基因**中有多少存在多个ployA和TSS的

> 大约有94.9%和98.7%的多转录本基因，存在多个TSS和ployA位点
>
> 平均每个基因存在TSS和ployA的数目：
>
> * A2   2.8个TSS和3.0个ployA
> * D5   2.4个TSS和2.7个ployA
> * TM1 2.27个TSS和2.4个ployA

| 基因组 | 多个TSS       | 多个ployA基因数  | 单个TSS | 单个ployA |
| --- | ----------- | ----------- | ----- | ------- |
| A2  | 13447 94.9% | 13994 98.7% | 720   | 173     |
| D5  | 11195 93.9% | 11710 98.3% | 716   | 201     |
| TM1 | 17448 93.7% | 18289 98.3% | 1158  | 317     |

**举一个重要基因,存在多个ployA和TSS的例子**，设计RACE进行扩增。

> TSS和ployA产生不同的转录本：两个转录本之间剪接位点完全相同，只在ployA或者TSS位点上存在差异的转录本数。**由于ployA和TSS的不同造成转录本间差异的比例**。

| 基因组 | 剪接位点都相同但TSS或者ployA位点不同 |
| --- | ---------------------- |
| A2  | 19794(29.5%)           |
| D5  | 15242(29.3%)           |
| TM1 | 20189(24.2%)           |

## **分析测序深度对于基因转录的影响**

测序饱和度分析,统计每个转录本的FL-read数目

统计每个转录本对应的基因数

```bash
##将注释文件中所有的转录本和基因信息提取出来，包括scaffold的，剔除比对到非基因区域的转录本
awk  '$3~/transcript/{print $0}' merge.gtf |awk  -F ";" '{print $3"\t"$2}'|sed -e 's/orginal_gene_id //g' -e 's/transcript_id //g' -e 's/\"//g' |sed 's/ //g' |awk '$1!~/PB/&&$2~/PB/{print $0}' >../../02Saturation_cruve/gene_Isoform.txt
##检测到参考基因注释的基因数和转录本对应的fl-read数
python ~/github/zpliuCode/Isoseq3/01saturation_curve/stattic_PacBio_readCount.py  -quiver  ../CO31_32_result/04_Cluster/quivered.all.fa   -collapse ../CO31_32_result/06_Alignment/all.collapsed.group.txt -gene ./gene_Isoform.txt  -out gene_isoform_readCount.txt
```

进行饱和曲线分析

* 根据测得的FL-read数目进行分析
* 依次随机10%、20%\~100%的FL-read；看会有多少基因和转录本被测到

```bash
##统计FL read 数目
grep ">" ../CO31_32_result/03_Classify/total.flnc.fasta |wc -l
##随机取样
 python ~/github/zpliuCode/Isoseq3/01saturation_curve/saturation_cruve.py  TM1_gene_isoform_readCount.txt  479445 plot_gene_isoform_saturation_cruve.txt
 awk '{print $0"\tfullMatch"}' full_match_cruve.txt  >>TM1_plot_cruve.txt
 awk '{print $0"\tincompleteMatch"}' incomplete_spliceMatch_cruve.txt >>TM1_plot_cruve.txt
 awk '{print $0"\tunAnnotion"}' unAnnotion_cruve.txt >>TM1_plot_cruve.txt
```

转录本长度和支持的FL-read 数目热图。