设计引物
之前看用RNA-seq去看转录本的表达量的的时候,由于模板浓度太低了,跑不出发生AS的那条带;因此使用PacBio的read数目来进一步筛选发生剪接的转录本。
#A2 vs At
evolution4/A2_vs_AT/Isoform_DomainFPKMChange/IsoformDomain_readCount.txt
##查看read数目和基因数目
grep -v "Nochange" IsoformDomain_readCount.txt |cut -f1,3,8,9|head
#D5 vs Dt
使用发生AS的基因和productive transcript发改变的基因取交集,同时满足PacBio reads数目大于10
grep Nochange -v conserve_specialIsoformDomain.txt |cut -f3|cat - AS_gene.txt |sort |uniq -d|xargs -I {} grep {} ../../../geneFPKM_readsCount/TM1_geneFPKM_read.txt
比较不同棉种中productive转录本序列是否发生较大的缺失
当提取两个棉种中主要表达的转录本进行,预测它们的CDS的时候;有时候它们的cDNA仅仅只相差个几bp,都可能会影响后续CDS的预测;因此还要进行进一步筛选,去比较两个cDNA是否存在较大的indel
脚本用法: https://github.com/BiocottonHub/zpliuCode/tree/master/transcriptSV#transcriptsv
##鉴定祖先和后代中主要表达的转录本是否存在PAV
python ~/github/zpliuCode/transcriptSV/transcriptSvs.py -homolog ~/work/Alternative/result/homologo/homologGene/Result/A2_vs_At_collinearity.txt -fasta1 ~/work/Alternative/result/Ga_result/CO11_12_result/collapse/mappingTogene.fa -fasta2 ../../collapse/mappingTogene.fa -RNAseq1 ../A2_PacBio.txt -RNAseq2 ../TM1_PacBio.txt -out ./A2_vs_At_productive_isoform.txt
是否存在AS导致的productive 转录本发生改变。
这里主要表达的转录本还是看full-read 数目是不是最多的。
首先提取同源基因中主要表达的转录本
比较两者是否存在缺失。
比较两者CDS序列是否存在差异
存在以下几种情况:
两个转录本CDS序列一致,而转录本存在结构变异;也就是AS影响的是OUTR区域。
两个转录本CDS序列不一致,并且转录本存在结构变异;说明可能是AS影响了编码区
两个转录本的CDS序列一致,同时转录本不存在结构变异;说明productive 转录本是保守的。
两种转录本的CDS序列不一致,而不存在结构变异,这就可能由于MUMmer只能鉴定比较长的PAV
例如Gorai.010G219800与Ghir_D06G020630相比存在一个主要表达转录本的差异,在D5、At中第一个内含子都是组成型的剪接,而在Dt中第一个外显子持续性的保持外显子的状态;这还是个热激蛋白。
从这个例子中,我们不难发现有的时候D5、Dt可能只转录出了一种类型的转录本,因此没能检测到可变剪接的存在,但将D5与Dt转录出的转录本进行比较,我们可以发现它们的转录本类似的发生了AS。
TM-1
A2
D5
Last updated
Was this helpful?