全长转录本数据的统计

全长转录本数据统计

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

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

 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

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，说明基因在转录过程中存在丰富的剪接形式。

与此同时两个二倍体物种相比于四倍体物种有着更加多的转录本

##统计每个基因的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

##提取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组成的。

• 单外显子转录本

• 多外显子基因，其中多少基因会产生多个转录本结构

• 分析多外显子转录本的长度分布

• 分析参考基因组中鉴定到的多个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

##提取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

##比较基因内转录本间的剪接位点和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数目

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

##将注释文件中所有的转录本和基因信息提取出来，包括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；看会有多少基因和转录本被测到

##统计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 数目热图。