根据bulk RNA-seq或者single-cell RNA-seq来找isoform的类型。

 

难点：二代NGS的reads长度较短，可能无法成功组装出full-length的isoform。

 

我们只需要关注junction read，基数即可。

 

junction read：在成熟mRNA中，测序的reads如果同时跨越了两个exon，则为junction read。

如何提取出junction read？

 

/home/lizhixin/softwares/regtools/build/regtools junctions extract -r chr19:797075-812327 -s 0 UE.bam

　　

自己写一个程序吧，用bedtools coverage这个功能。

bedtools coverage -a regions.bed -b reads.bam

　　

spliceGraph

 

核心就是count reads，同时比对到两个点上。【这个不难，但是并不能说明问题】

samtools view -u -@ 2 UE.bam chr19:803636-803636 | less -S

　　

 必须鉴定出拼接两个exon的read，这样的才是有意义的证据，单纯的连接的数据并不实用【中间可能插入了很多个exon，连接≠拼接】。

 

cat gencode.v27.annotation.gtf | grep "\"PKM\"" | awk '{if($3=="exon")print $0}' > PKM.txt

　　

 

samtools view -b -@ 2 23c9.bam chr19:797075-812327 > PTBP1.23c9.bam

　　

 

参考：

Question: Count Junctions In A Sam/Bam File

junctions extract - RegTools

Question: How Many Reads In A Bam File Are Aligned To a Specific Region