重新下载原始数据进行比对

• Ga SRR5886156

• Gr SRR5886149

• Gh SRR5886147

hisat2比对

• -novel-splicesite-outfile生成剪切位点文件

• --max-intronlen指定最大内含子长度

• --rf建库方式

##构建索引

bsub -J TM1 -n 10 -o %J.hisat2.out -e %J.hisat2.err -R span[hosts=1]  "hisat2 -x ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/hisat2_build/TM1.0 -1 ../Trimmomatic/SRR5886147_Trimm_1P  -2 ../Trimmomatic/SRR5886147_Trimm_2P -p 10 --novel-splicesite-outfile TM1_splicesite.txt  --max-intronlen 10000 --rf  -S TM1_leaf.sam"

## 过滤和去重PCR重复
bsub -J samtools -n 2 -o %J.samtools.out -e %J.samtools.err -R span[hosts=1] "samtools view -q 20 -S TM1_leaf.sam  -@ 2 -b -o TM1_leaf.bam && samtools sort -@ 2 TM1_leaf.bam -O bam -o TM1_sorted.bam &&samtools rmdup TM1_sorted.bam TM1_rmdup.bam"

Stringtie计算gene表达量

stringtie hisat2/TM1_rmdup.bam -G ~/work/Alternative/result/homologo/FEST3/isforms/TM1/geneExpress/isform.gtf  -b ./stringtie/  -p 10

计算PIR、PSI

PIR：IR的比例

PSI：intron被剪切的比例

将RNA-seq比对文件转为bed

BAM文件中read比对说明

M 完全匹配

I 插入

D 缺失

N gap

S 替代

• -cigar输出每个read的比对情况

• -split不会连接有gap的read

80M20N30M这种评分则会将这个read输出两个bed坐标

~/software/bedtools2-2.29.0/bin/bamToBed -cigar -split -i TM1_rmdup.bam >TM1.bed
## 过滤掉scaffold上的结果
awk '$1~/^Ch/{print $0}' D5.bed  >11 
##提交任务
bsub -q q2680v2 -n 20 -J A2 -e %J.err -o %J.out -R span[hosts=1] "python ~/scripte/Alternative/module/FEST3/PSI/IR_PSI.py  -AS ~/work/Alternative/result/Ga_result/CO11_12_result/11_AS/end_splice.txt3 -bed ../geneExpress/hisat2/A2.bed  -o A2_PSI.txt -p 20 "

计算所有Intron的PSI

• 合并所有的intron

• 计算每个intron的PIR值

## 提取intron坐标
python extractAllIntron.py  -pgtf ~/work/Alternative/result/Gr_result/CO41_42_result/07_annotation/merge.gtf -rgtf ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gtf -o 11
awk '$1~/^C/{print $1"\t"$2"\t"$3}' 11 |sort -k1,1 -k2,3n|uniq  >D5_intron.txt
## 合并后计算每个intron的PIR值

AS剪切在一定程度上造成了isform的多样性

发生剪切的高频区域，每个转录本所包含AS 事件的数目

比较由AS产生的保守isform发生NMD的比例、和由AS产生的uniq isform发生NMD的比例

造成uniq isform表达水平比conserve isform表达水平低的原因，可能是NMD途径导致uniq isform降解；

先把所有转录本的ORF预测了

• 需要反转模式预测reverse

/public/home/software/opt/bio/software/EMBOSS/6.5.7/bin/getorf   -find 3 -noreverse -sequence cDNA.fa

根据终止密码子与exon-exon junction的位置判断是否会发生NMD

在处理D5时，有部分基因没有exon注释，搞出了，修改gtf文件然后再跑一次，还有D5参考基因组内没有cDNA序列，需要提取一次

我对读取gtf文件的类做了补充，支持没有exon注释的ISOfrom

• 修改D5 gtf文件

• 提取D5 cDNA序列

##提取参考基因组中cDNA序列
python ~/scripte/Alternative/extractcDNA.py -gtf ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gtf -fasta ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.0.fa -cDNA D5_cDNA.fa

#判断isoform是否NMD
python ~/scripte/Alternative/module/FEST3/NMD.py  -gtf ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gtf -orf 010g215500.orf  -o 11  >222
##统计多少是NMD的比例
awk '$4==1{a+=1}END{print a/NR}' ref_NMD.txt

关于对参考基因组中注释的isform的验证，需要筛选一下

• 筛选那些每个基因中表达量大于100 fpkm的注释的转录本

• 使用参考基因组的注释信息去跑一下这些转录本的表达量

we screened a set of 4,868 Arabidopsis genes, for which the annotated longest ORF is fully supported by cDNA evidence

cDNA序列序列5’和3‘已经确定了，感觉不需要反转

对鉴定到的PacBio isform进行NMD分析

• 距离 最后一个exon-exon的junction大于50bp就认为是潜在的NMD

## 对参考基因组的转录本和PacBio注释的转录本，预测NMD
python ~/scripte/Alternative/module/FEST3/NMD.py  -gtf merge.gtf  -orf merge.orf -o PacBio_NMD.txt

棉种	PacBio isform	NMD
TM1	89411	35218
A2	72374	26583
D5	55365	19330

N5xSMV.png

比较可能发生NMD的isform与没有NMD的表达水平差异

##多个棉种画同一个图里比较FPKM的话，量化一下；以表达量最多的那个搞
TM1统一提高1.35倍
##所有PacBio isoform的表达水平
awk 'NR>1{print $6"\t"$12}' ../../../geneExpress/stringtie/A2/t_data.ctab   >PacBio_FPKM.txt
##提取NMD isform的表达水平
grep "PB" PacBio_NMD.txt |awk '$4==0{print $1}'|cat - PacBio_FPKM.txt |sort -k1,1  -k2,2nr|awk '$2==""{print "0\t"$1}$2!=""{print $2"\t"$1}'|uniq  -f1 -d|awk '{print $1"\tD5\tNMD"}' >>compare_FPKM.txt
##提取none-NMD isform表达水平
 grep "PB" PacBio_NMD.txt |awk '$4!=0{print $1}'|cat - PacBio_FPKM.txt |sort -k1,1  -k2,2nr|awk '$2==""{print "0\t"$1}$2!=""{print $2"\t"$1}'|uniq  -f1 -d|awk '{print $1"\tD5\tnone"}' >>compare_FPKM.txt

分析AS与NMD

从剪切文件中提取两个转录本之间只有一种AS的差异，对没有发生AS的那个转录本进行ORF的预测，同时也对发生了AS的isoform预测ORF，选取相同的ST位点，看没有发生AS的转录本是否存在PTC并且PTC与exon-exon junction大于50bp；分析每种剪切类型中存在NMD的转录本

• 获取两个转录本只存在一种类型AS的差异

• 选择没有发生AS的转录本预测的ORF作为参考

• 计算参考转录本所预测的ST和TT

• 计算发生与参考转录本ST相同的ORF，比较TT是否提前

##合并两次预测的orf
grep ">P" 1.orf >merge.orf
grep ">G" ghir_a01g000010.orf >>merge.orf 
##合并两个gtf文件
cat ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf ../../../isforms/TM1/geneExpress/isform.gtf >merge.gtf
##分析
 python ~/scripte/Alternative/module/FEST3/NMD2.py -gtf merge.gtf  -orf merge.orf  -AS ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt3  -o 1111111
##判断是否有PTC
grep "ExonS" 1111111 |awk '$4<$5&&$5<$7{print $0}$4>$5&&$5>$7{print $0}'|wc -l
##对每种类型计算对应的NMD数目和比例
grep "IntronR" AS_NMD.txt |awk '$4<$5&&$5<$7{a+=1}$4>$5&&$5>$7{a+=1}END{print a"\t"NR"\t"a/NR}'

基因组	IR	ES	AltA	AltD
TM1	15100/18733	758/1274	1960/2825	1529/2216
A2	12887/15287	499/727	1262/1790	1029/1406
D5	10051/12296	675/917	1342/1843	980/1383

不同AS产生的NMD，在表达水平上的差异

##提取不同AS产生的isoform id
for i in IntronR ExonS AltA AltD ; do grep $i AS_NMD.txt |awk '$4<$5&&$5<$7{print $2}$4>$5&&$5>$7{print $2}'|cat - PacBio_FPKM.txt |sort -k1,1  -k2,2nr|awk '$2==""{print "0\t"$1}$2!=""{print $2"\t"$1}'|uniq  -f1 -d|awk '{print $1"\tD5\t""'$i'"}' >>As_FPKM.txt; done

N5xJzt.png