可变剪切 | isoform | 提取特定exon的usage

可变剪切 | isoform | 提取特定exon的usage | DEXSeq

参考：Inferring differential exon usage in RNA-Seq data with the DEXSeq package

http://bioconductor.org/packages/release/bioc/html/DEXSeq.html 【教程已经非常清楚，需要耐心仔细的研读】

跟着tutorial走一遍【生信码农的基本功：工具测试评比】

关键问题：

gtf文件 - 必须是来自Ensembl的，否则没法正确运行，版本和格式是永恒的浪费时间的点
单细胞bam - 必须合并，不然出来的都是0，全部都检测不到
运行速度 - Python的wrapper是真的慢

更新：

之前一直模糊的概念，这里必须搞清楚！！！

read count

每个exon or exonic part的read count
每个样本的read count总和，即library size 测序深度
每个gene的read count总和，即该基因的所有transcript的read总数
relative usage of an exon，即 this / others，该exon的read总数除以比对到该gene其他exon的read总数 = 相对exon usage

问题：

把count数据转成CPM就能比较exon的CPM吗？ - 不能，exon层面就不能简单用gene那一套了，必须先关注这个gene本身是否有差异表达
relative usage of an exon的数据可比吗？ - 这就是校正了gene表达总量后的值，用来比较就更有意义了。

参考：

关于normalization - Introduction to DGE - ARCHIVED

Gene length corrected trimmed mean of M-values (GeTMM) processing of RNA-seq data performs similarly in intersample analyses while improving intrasample comparisons

先拆

ls archive_20201222/bam/BAM.sortedByCoord.out/*.bam | sed "s:^:`pwd`/: "  > all.bam.csv

cat all.bam.csv | grep "10c2" > 10c2.list
cat all.bam.csv | grep "17C8" > 17C8.list
cat all.bam.csv | grep "1C11" > 1C11.list
cat all.bam.csv | grep "20c7" > 20c7.list
cat all.bam.csv | grep "23C9" > 23C9.list
cat all.bam.csv | grep "3C15" > 3C15.list
cat all.bam.csv | grep "5c3" > 5c3.list
cat all.bam.csv | grep "6c5" > 6c5.list
cat all.bam.csv | grep "/IMR90" > IMR90.list
cat all.bam.csv | grep "IMR-N_Diff-D20" > IMR-N_Diff-D20.list
cat all.bam.csv | grep "IMR-N_Diff-D40" > IMR-N_Diff-D40.list
cat all.bam.csv | grep "iPSC-IMR90" > iPSC-IMR90.list
cat all.bam.csv | grep "iPSC-UE" > iPSC-UE.list
cat all.bam.csv | grep "/UE" | grep -v "N_Diff" > UE.list
cat all.bam.csv | grep "UE-N_Diff_D20" > UE-N_Diff_D20.list
cat all.bam.csv | grep "UE-N_Diff-D40" > UE-N_Diff-D40.list

split -d -n l/5 10c2.list 10c2.list
split -d -n l/5 17C8.list 17C8.list
split -d -n l/5 1C11.list 1C11.list
split -d -n l/5 20c7.list 20c7.list
split -d -n l/5 23C9.list 23C9.list
split -d -n l/5 3C15.list 3C15.list
split -d -n l/5 5c3.list 5c3.list
split -d -n l/5 6c5.list 6c5.list
split -d -n l/5 IMR90.list IMR90.list
split -d -n l/5 IMR-N_Diff-D20.list IMR-N_Diff-D20.list
split -d -n l/5 IMR-N_Diff-D40.list IMR-N_Diff-D40.list
split -d -n l/5 iPSC-IMR90.list iPSC-IMR90.list
split -d -n l/5 iPSC-UE.list iPSC-UE.list
split -d -n l/5 UE.list UE.list
split -d -n l/5 UE-N_Diff_D20.list UE-N_Diff_D20.list
split -d -n l/5 UE-N_Diff-D40.list UE-N_Diff-D40.list

合并

samtools merge -b UE.bam.list -@ 2 -r -p UE.bam
samtools merge -b UE.bam.list -@ 2 -r -p UE.bam

批量【先合并两个测试一下】

for i in `ls {IMR90.list0*,17C8.list0*}`
do
echo $i
samtools merge -b $i -@ 2 -r -p ${i}.bam &&
echo "done!!!"
done

准备特殊的注释文件 & 计数

gencode版本

# gencode
# ~/databases/hg19/
python ~/softwares/R_lib_361/DEXSeq/python_scripts/dexseq_prepare_annotation.py gencode.v27.annotation.gtf gencode.v27.annotation.DEXSeq.chr.gff

# ~/project/scRNA-seq/rawData/smart-seq/archive_20201222/bam/BAM.sortedByCoord.out/IMR90-E2-A17_STAR_Aligned.sortedByCoord.out.bam
python ~/softwares/R_lib_361/DEXSeq/python_scripts/dexseq_count.py -p yes -f bam ~/databases/hg19/gencode.v27.annotation.DEXSeq.chr.gff ~/project/scRNA-seq/rawData/smart-seq/merged.bams/bam/IMR90.list00.bam tmp.result.txt

ensemble版本

############################################
# ensemble
# ~/references/SUPPA2/ref/hg19/
python ~/softwares/R_lib_361/DEXSeq/python_scripts/dexseq_prepare_annotation.py Homo_sapiens.GRCh37.75.formatted.gtf Homo_sapiens.GRCh37.75.DEXSeq.chr.gff

python ~/softwares/R_lib_361/DEXSeq/python_scripts/dexseq_count.py -f bam ~/references/SUPPA2/ref/hg19/Homo_sapiens.GRCh37.75.DEXSeq.chr.gff ~/project/scRNA-seq/rawData/smart-seq/merged.bams/bam/IMR90.list00.bam tmp.result.txt

批量运行

# 批量运行
for i in `ls ~/project/scRNA-seq/rawData/smart-seq/merged.bams/bam/*.bam`
do
python ~/softwares/R_lib_361/DEXSeq/python_scripts/dexseq_count.py -f bam ~/references/SUPPA2/ref/hg19/Homo_sapiens.GRCh37.75.DEXSeq.chr.gff $i ${i}.result.txt &&
echo $i
done

文件如下：

chr15   dexseq_prepare_annotation.py    aggregate_gene  72491370        72524164        .       -       .
       gene_id "ENSG00000067225"
chr15   dexseq_prepare_annotation.py    exonic_part     72491370        72491372        .       -       .
       gene_id "ENSG00000067225"; transcripts "ENST00000564440+ENST00000319622"; exonic_part_number "001"
chr15   dexseq_prepare_annotation.py    exonic_part     72491373        72491376        .       -       .
       gene_id "ENSG00000067225"; transcripts "ENST00000319622+ENST00000565184+ENST00000389093+ENST00000335181+ENST00000568883+ENST00000564440"; exonic_part_number "002"
chr15   dexseq_prepare_annotation.py    exonic_part     72491377        72491445        .       -       .
       gene_id "ENSG00000067225"; transcripts "ENST00000319622+ENST00000565184+ENST00000389093+ENST00000335181+ENST00000565143+ENST00000568883+ENST00000564440"; exonic_part_number "003"
chr15   dexseq_prepare_annotation.py    exonic_part     72491446        72491930        .       -       .
       gene_id "ENSG00000067225"; transcripts "ENST00000567118+ENST00000319622+ENST00000565184+ENST00000389093+ENST00000335181+ENST00000449901+ENST00000565143+ENST00000568883+ENST00000564440"; exonic_part_number "004"

很明确，每一个都是pseudo exon【有时候会把exon切得非常碎，需要借助可视化工具来操作】

# check
# https://www.ncbi.nlm.nih.gov/gene/5315
PKM ENSG00000067225

总结：

对于大型bam文件来说，这个计数实在是太慢了，需要耐心等待
操作明确，对于pseudo exon进行read的计数，然后用已有的成熟工具来做差异分析
该工具画图居然需要root权限，实在是匪夷所思

可以作为一个AS分析的备选工具。

参考分析目录：

CGS: project/scPipeline/AS/DEXSeq/DEXSeq.ipynb

CGS: project/scRNA-seq/rawData/smart-seq/merged.bams/

待续~