ixxmu
commented
3 years ago 上游比对流程
Closed ixxmu closed 3 years ago
ixxmu
commented
3 years ago 上游比对流程
github-actions[bot]
commented
3 years ago 连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！大家开始根据我的ngs组学视频进行一系列公共数据集分析实战，其中一个小伙伴让我非常惊喜，不需要怎么沟通和指导，就默默的完成了一个实战！
下面是温州医科大硕士“何物昂”的笔记
数据集为GSE149638， 2x101 bp paired-end RNA-seq，Illumina HiSeq 2500 with poly-A selection。源于健康人的M0和M1 macrophages。原始数据M0和M1各有48个重复。全部使用还是需要耗费一定时间和计算资源的，这里就各挑选3个重复进行练习。
我比较喜欢去ebi里下载数据，感觉ebi下载数据更人性化一点。去ebi 搜索GSE149638
进入数据页面后，往下滑，是会看到read files的表格的，但根据个人网络，可能加载出来的比较慢。
因为当前默认显示的，不包含sample_title，其余信息好像也没法直接区分M0和M1数据。我们可以借助sample title 来区分M0和M1数据。fastq_ftp 的url可以放到迅雷里去下载，fastq_aspera可以用aspera下载。这次练习我就用aspera下载。
下载后，会得到filereport_read_run_PRJNA629498_tsv.txt 这个文件，里面的内容是这样子的。
如果是用迅雷下载数据的话，参考下面这命令修改下。
grep "HMDM M1" filereport_read_run_PRJNA629498_tsv.txt | cut -f 7 | head -n 3 | sed -r 's/;/\n/' | awk '{print "ftp://"$0}'
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/032/SRR11652532/SRR11652532_1.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/032/SRR11652532/SRR11652532_2.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/034/SRR11652534/SRR11652534_1.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/034/SRR11652534/SRR11652534_2.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/036/SRR11652536/SRR11652536_1.fastq.gz
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR116/036/SRR11652536/SRR11652536_2.fastq.gz把上面输出的url复制到迅雷里就可以下载了
命令解释：
grep "HMDM M1" 用于查询M1数据
cut -f 7 取第7列，当前我的第7列是ftp url，不过这并不固定，取决于之前勾选的columns。
head -n 3 取前3列，即M1的3个重复
sed -r 's/;/\n/' PE 数据r1和r2数据，是以";"隔开，这里将";"替换成换行符,SE不需要
awk '{print "ftp://"$0}' 添加ftp://头，迅雷才可以直接识别下载
使用aspera下载 ，具体参考 使用aspera从EBI下载fastq数据，抛弃NCBI的SRA数据库
python dl_fq.py filereport_read_run_PRJNA629498_tsv.txt -n 6 --id_rsa ~/.aspera/connect/etc/asperaweb_id_dsa.openssh
filereport_read_run_PRJNA629498_tsv.txt 是从ebi里下载得来的，
-n 是指定下载多少个样本（不区分M0和M1的）
--id_rsa ~/.aspera/connect/etc/asperaweb_id_dsa.openssh 指定密钥位置
运行过程apsera 输出日志
下载完成查看，名字已经改了。
ls -lh |cut -d" " -f 5-
2.3G Jun 25 13:55 HMDM_M0_rep1_SRR11652531_1.fastq.gz
2.2G Jun 25 14:04 HMDM_M0_rep1_SRR11652531_2.fastq.gz
2.1G Jun 25 14:34 HMDM_M0_rep2_SRR11652533_1.fastq.gz
2.1G Jun 25 14:41 HMDM_M0_rep2_SRR11652533_2.fastq.gz
2.7G Jun 25 15:04 HMDM_M0_rep3_SRR11652535_1.fastq.gz
2.7G Jun 25 15:16 HMDM_M0_rep3_SRR11652535_2.fastq.gz
2.3G Jun 25 14:14 HMDM_M1_rep1_SRR11652532_1.fastq.gz
2.3G Jun 25 14:23 HMDM_M1_rep1_SRR11652532_2.fastq.gz
2.7G Jun 25 14:48 HMDM_M1_rep2_SRR11652534_1.fastq.gz
2.6G Jun 25 14:55 HMDM_M1_rep2_SRR11652534_2.fastq.gz
2.5G Jun 25 16:20 HMDM_M1_rep3_SRR11652536_1.fastq.gz
2.5G Jun 25 16:00 HMDM_M1_rep3_SRR11652536_2.fastq.gz
运行过程中，可能由于网络原因，下载失败，再次运行命令就好了，已经下载好的，不会重复下载的。
首先，做质控分析，看下fq文件质量情况. 使用软件为fastqc,multiqc. fastqc用于做质控的，multiqc用于整合多个样本质控结果的。
安装都可以通过conda
conda install -c bioconda fastqc -y
conda install -c bioconda multiqc -y
运行命令
mkdir fastqc
fastqc *fastq.gz -o fastqc/ -t 10
multiqc fastqc/ -o multiqc -n HMDM
fastqc 参数：
*fastq.gz 匹配当前目录下的以fastq.gz结尾的fq文件
-o fastqc 指定输出为 fastqc，需提前创建
-t 10 指定线程数
multiqc 参数
fastqc/ : fastqc 程序的输出目录
-o multiqc : multiqc 程序的输出目录
-n HMDM ：指定输出文件前缀名
查看程序输出
ll fastqc/ | cut -d " " -f 5-
240544 Jun 25 17:23 HMDM_M0_rep1_SRR11652531_1_fastqc.html
268242 Jun 25 17:23 HMDM_M0_rep1_SRR11652531_1_fastqc.zip
237383 Jun 25 17:22 HMDM_M0_rep1_SRR11652531_2_fastqc.html
266669 Jun 25 17:22 HMDM_M0_rep1_SRR11652531_2_fastqc.zip
239187 Jun 25 17:22 HMDM_M0_rep2_SRR11652533_1_fastqc.html
269215 Jun 25 17:22 HMDM_M0_rep2_SRR11652533_1_fastqc.zip
......
257073 Jun 25 17:24 HMDM_M1_rep2_SRR11652534_2_fastqc.zip
237224 Jun 25 17:23 HMDM_M1_rep3_SRR11652536_1_fastqc.html
263355 Jun 25 17:23 HMDM_M1_rep3_SRR11652536_1_fastqc.zip
234936 Jun 25 17:23 HMDM_M1_rep3_SRR11652536_2_fastqc.html
259976 Jun 25 17:23 HMDM_M1_rep3_SRR11652536_2_fastqc.zip
ll multiqc/ | cut -d " " -f 5-
160 Jun 25 19:37 HMDM_data
1294537 Jun 25 19:37 HMDM.html
可以看到fastqc程序是每个fq文件对应一个质量报告，使用multiqc可以将他们整合在一起。用浏览器打开html文件就可以查看分析报告了。
下面进行接头修剪（虽然我查看了，上面的6个样本，没检测出接头,不过还是跑下流程）
adapter 去除使用trim_galore， 该软件可以自动识别接头。
安装
conda install -c bioconda trim-galore
程序运行
mkdir ../cleanData
ls *gz|cut -d"_" -f 4 | sort -u | \
while read id;do \
nohup trim_galore -q 25 --phred33 --length 36 \
--stringency 3 --paired -o ../cleanData *${id}*.gz & \
done
因为这里样本量少，就直接 nohup到后台了。
输出结果：
ll ../cleanData/*fq.gz | cut -d " " -f 5-
2361440662 Jun 25 18:08 ../cleanData/HMDM_M0_rep1_SRR11652531_1_val_1.fq.gz
2299024844 Jun 25 18:08 ../cleanData/HMDM_M0_rep1_SRR11652531_2_val_2.fq.gz
2109567453 Jun 25 18:02 ../cleanData/HMDM_M0_rep2_SRR11652533_1_val_1.fq.gz
2097839448 Jun 25 18:02 ../cleanData/HMDM_M0_rep2_SRR11652533_2_val_2.fq.gz
2759870348 Jun 25 18:19 ../cleanData/HMDM_M0_rep3_SRR11652535_1_val_1.fq.gz
2785232279 Jun 25 18:19 ../cleanData/HMDM_M0_rep3_SRR11652535_2_val_2.fq.gz
2345725234 Jun 25 18:09 ../cleanData/HMDM_M1_rep1_SRR11652532_1_val_1.fq.gz
2392756954 Jun 25 18:09 ../cleanData/HMDM_M1_rep1_SRR11652532_2_val_2.fq.gz
2737439459 Jun 25 18:17 ../cleanData/HMDM_M1_rep2_SRR11652534_1_val_1.fq.gz
2698516298 Jun 25 18:17 ../cleanData/HMDM_M1_rep2_SRR11652534_2_val_2.fq.gz
2536073207 Jun 25 18:12 ../cleanData/HMDM_M1_rep3_SRR11652536_1_val_1.fq.gz
2524396378 Jun 25 18:12 ../cleanData/HMDM_M1_rep3_SRR11652536_2_val_2.fq.gz
修剪完，再做质控。
cd ../cleanData
mkdir fastqc
fastqc *fq.gz -o fastqc/ -t 10
multiqc fastqc/ -o multiqc -n HMDM
这里采用salmon来进行转录本定量，比传统的使用STAR/HISAT2等比对软件比对参考基因组，然后进行raw count计数更快。使用salmon也需要先构建index。salmon提供了三种构建方式。
salmon可以只用转录组序列来构建index,这样也更快。不过作者推荐使用全基因组数据，这可以得到更好的效果。这里我们就采用全基因组来构建索引。从这里下载基因组，转录组。
http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/
转录组: http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.transcripts.fa.gz
基因组：http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/GRCh38.primary_assembly.genome.fa.gz
salmon 安装
conda install -c bioconda salmon
salmon index 构建
cd ../ref # 序列已经下载到ref目录下了
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) | cut -d " " -f 1 > decoys.txt
sed -i.bak -e 's/>//g' decoys.txt
cat gencode.v38.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz
salmon index -t gentrome.fa.gz -d decoys.txt -p 12 -i salmon_index --gencode
这步会花比较长的时间。闲慢的话，考虑直接使用转录组进行index
salmon index -t gencode.v38.transcripts.fa.gz -i salmon_index
salmon 定量
cd ../cleanData
mkdir ../salmon_output
index=../ref/salmon_index
ls *gz|cut -d "_" -f 1,2,3,4 | sort -u | while read id;do \
nohup salmon quant -i $index -l A --gcBias \
-1 ${id}*_1_val_1.fq.gz -2 ${id}*_2_val_2.fq.gz -p 2 \
-o ../salmon_output/${id}_output 1>${id}_salmon.log 2>&1 & done
salmon参数说明
-i $index : $index 是 ../ref/salmon_index 之前构建的索引
-1， -2 ：指PE 的r1,r2
-p 2: 指定线程数
-o ：输出目录
--gcBias：用于自动矫正GC偏好性
-l A ：--libType 指测序文库类型，分单双端，是否链特异，read方向。salmon指定了一些字符来描述这些不同的文库类型。参数值A 指让salmon自动推测文库类型。
查看输出文件,其中主要用的信息在quant.sf文件里
ll ../salmon_output/HMDM_*_output/quant.sf | cut -d " " -f 5-
10912133 Jun 25 21:08 ../salmon_output/HMDM_M0_rep1_SRR11652531_output/quant.sf
10899997 Jun 25 21:00 ../salmon_output/HMDM_M0_rep2_SRR11652533_output/quant.sf
10924786 Jun 25 21:17 ../salmon_output/HMDM_M0_rep3_SRR11652535_output/quant.sf
10910812 Jun 25 21:11 ../salmon_output/HMDM_M1_rep1_SRR11652532_output/quant.sf
10913099 Jun 25 21:11 ../salmon_output/HMDM_M1_rep2_SRR11652534_output/quant.sf
10914575 Jun 25 21:10 ../salmon_output/HMDM_M1_rep3_SRR11652536_output/quant.sf
head ../salmon_output/HMDM_M0_rep1_SRR11652531_output/quant.sf
Name Length EffectiveLength TPM NumReads
ENST00000456328.2 1657 1479.000 0.000000 0.000
ENST00000450305.2 632 454.000 0.000000 0.000
ENST00000488147.1 1351 1116.200 3.598279 126.358
ENST00000619216.1 68 5.000 0.000000 0.000
ENST00000473358.1 712 534.000 0.000000 0.000
ENST00000469289.1 535 358.000 0.000000 0.000
ENST00000607096.1 138 17.000 0.000000 0.000
ENST00000417324.1 1187 1009.000 0.000000 0.000
ENST00000461467.1 590 412.000 0.000000 0.000
salmon定量出来的结果不能直接用于差异基因的寻找，需要转换成基因层面count matrices。我们使用tximport来进行转换。R里写代码处理，使用DESeq2进行差异基因寻找。
library(GenomicFeatures)
library(AnnotationDbi)
library(tximport)
library(DESeq2)
### 我们直接从gtf注释文件（和基因组转录组对应版本的注释）获得基因和转录本的对应关系
hg38_txdb <- makeTxDbFromGFF("../ref/gencode.v38.annotation.gtf",
format="gtf",
organism="Homo sapiens",
dbxrefTag = "gene_name")
k <- keys(hg38_txdb, keytype = "TXNAME")
tx2gene <- select(hg38_txdb, k, "GENEID", "TXNAME")
head(tx2gene)
# TXNAME GENEID
#1 ENST00000456328.2 ENSG00000223972.5
#2 ENST00000450305.2 ENSG00000223972.5
#3 ENST00000473358.1 ENSG00000243485.5
#4 ENST00000469289.1 ENSG00000243485.5
#5 ENST00000607096.1 ENSG00000284332.1
#6 ENST00000606857.1 ENSG00000268020.3
quant.files <- list.files("../salmon_output", pattern = "quant.sf",
full.names = T, recursive = T)
txi <- tximport(quant.files, type = "salmon", tx2gene = tx2gene)
names(txi)
# [1] "abundance" "counts" "length" "countsFromAbundance"
sampleTable <- data.frame(condition = factor(rep(c("M0", "M1"), each = 3)))
rownames(sampleTable) <- paste0("sample", 1:6)
### DESeq2提供了直接从Tximport构建数据的函数
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
## 过滤一下，6个样本总counts数大于60
keep <- rowSums(counts(dds)) >= 60
dds <- dds[keep,]
### pca 作图
library(ggfortify)
## normalization
count.sf <- estimateSizeFactors(dds)
normal.mtx <- counts(count.sf, normalized=T)
df_all <- cbind(t(normal.mtx), data.frame(group=c(rep('M0', 3),rep('M1', 3))))
pca_data <- prcomp(t(normal.mtx ), scale. = F)
autoplot(pca_data, data = df_all, colour = 'group', label = TRUE) + theme_bw()
pca 图
### 差异分析
deseq2.obj <- DESeq(dds)
# 获取结果
res <- results(deseq2.obj)
res.df <- as.data.frame(res)
r$> head(res.df)
baseMean log2FoldChange lfcSE stat pvalue padj
ENSG00000000003.15 12.78531 -1.15948975 1.0905031 -1.0632613 2.876635e-01 3.982777e-01
ENSG00000000419.14 1153.61638 -0.05190081 0.2331701 -0.2225878 8.238563e-01 8.760466e-01
ENSG00000000457.14 842.03593 2.47434389 0.3036221 8.1494188 3.656775e-16 8.626111e-15
ENSG00000000460.17 264.26436 -1.72342072 0.4529815 -3.8046160 1.420244e-04 5.627837e-04
ENSG00000000938.13 6142.70208 -2.29502017 0.3162089 -7.2579238 3.930770e-13 6.552265e-12
ENSG00000000971.16 1376.49691 4.16870903 0.5808924 7.1763881 7.157718e-13 1.156001e-11
r$> keep <- abs(res$log2FoldChange) > 1.5 & res$padj < 0.01
r$> deg <- res.df[keep, ]
r$> dim(deg)
[1] 4083 6
## 使用 https://mp.weixin.qq.com/s/NC5mKwAsYbm1Q-PFC5IEOg 代码
library(pheatmap)
dat <- normal.mtx
FC <- res.df$log2FoldChange
names(FC) <- rownames(res.df)
# 排序差异倍数，提取前100和后100的基因名
DEG_200 <- c(names(head(sort(FC),100)),names(tail(sort(FC),100)))
# 提取基因的归一化
dat <- t(scale(t(dat[DEG_200,])))
group_list <- df_all$group
# 添加注释条
group <- data.frame(group=group_list)
rownames(group) <- colnames(dat)
# 大于2的值赋值为2
dat[dat > 2] <- 2
# 低于-2的值赋值为-2
dat[dat < -2] <- -2
pheatmap(dat, cluster_cols = T,
show_colnames =F,show_rownames = F,
annotation_col=group,
color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
library(org.Hs.eg.db)
library(clusterProfiler)
## 获取上下调基因
deg.up <- res.df[res$log2FoldChange > 1.5 & res$padj < 0.01, ]
deg.down <- res.df[res$log2FoldChange < -1.5 & res$padj < 0.01, ]
## 去除ENSEMBL ID的版本号
deg.up.gene <- gsub("\\..*", "", rownames(deg.up))
deg.down.gene <- gsub("\\..*", "", rownames(deg.down))
deg.ls <- list(up = deg.up.gene, down = deg.down.gene)
ego.cmp <- compareCluster(deg.ls ,
fun = "enrichGO",
qvalueCutoff = 0.01,
pvalueCutoff = 0.05,
ont = 'BP',
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
)
dotplot(ego.cmp)
up.ego <- enrichGO(gene = deg.up.gene,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.01,
readable = TRUE,
pool = TRUE)
dotplot(up.ego, split="ONTOLOGY", showCategory = 15, title = "up") +
facet_grid(ONTOLOGY~., scale="free")
down.ego <- enrichGO(gene = deg.down.gene,
keyType = "ENSEMBL",
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.01,
readable = TRUE,
pool = TRUE)
dotplot(down.ego, split="ONTOLOGY", showCategory = 15, title = "down") +
facet_grid(ONTOLOGY~., scale="free")
### KEGG富集，需要使用ENTREZID
deg.up.gene.ei <- bitr(deg.up.gene, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
deg.down.gene.ei <- bitr(deg.down.gene, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
deg.ei.ls <- list(up = deg.up.gene.ei, down = deg.down.gene.ei)
kegg.cmp <- compareCluster(deg.ei.ls,
fun = "enrichKEGG",
qvalueCutoff = 0.01,
pvalueCutoff = 0.05,
organism = "hsa"
)
dotplot(kegg.cmp)
#### 单个富集查看的话，用函数enrichKEGG
可变剪切分析，这里使用SUPPA2 来进行。据作者说，这是一个又快，准确度又高，对样本测序深度要求还低的软件。
安装也可以通过conda
conda install -c bioconda suppa
PSI（proportion spliced-in）作为可变剪切中常见的评估指标，用于衡量可变剪切发生率。在不同软件中，计算可能有所不同。举一个skipping exon的例子，若F1 表示那些包含某个exon (这里研究local AS events, 指特定的一个exon)的转录本，F2 表示那些不包含那个exon 的转录本。则PSI,计算公式如下：
SUPPA2需要使用基因注释gtf文件中的exon注释信息来生成不同的可变剪切事件：
cd ../ref
wget http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.annotation.gtf.gz
gunzip gencode.v38.annotation.gtf.gz
suppa.py generateEvents -i gencode.v38.annotation.gtf -o suppa/gencode.v38.events -e SE SS MX RI FL -f ioe
ll -h suppa/ | cut -d " " -f 5-
17M Jun 26 13:34 gencode.v38.events_A3_strict.gtf
5.5M Jun 26 13:34 gencode.v38.events_A3_strict.ioe
15M Jun 26 13:34 gencode.v38.events_A5_strict.gtf
4.3M Jun 26 13:34 gencode.v38.events_A5_strict.ioe
103M Jun 26 13:34 gencode.v38.events_AF_strict.gtf
19M Jun 26 13:34 gencode.v38.events_AF_strict.ioe
34M Jun 26 13:34 gencode.v38.events_AL_strict.gtf
6.2M Jun 26 13:34 gencode.v38.events_AL_strict.ioe
14M Jun 26 13:34 gencode.v38.events_MX_strict.gtf
1.9M Jun 26 13:34 gencode.v38.events_MX_strict.ioe
5.0M Jun 26 13:34 gencode.v38.events_RI_strict.gtf
2.0M Jun 26 13:34 gencode.v38.events_RI_strict.ioe
49M Jun 26 13:34 gencode.v38.events_SE_strict.gtf
13M Jun 26 13:34 gencode.v38.events_SE_strict.ioe
# 若是想研究所有的AS事件类型，可以把所有的AS event ioe文件合并成一个
awk 'FNR==1 && NR!=1 { while (/^<header>/) getline; } 1 {print}' suppa/*.ioe > gencode.v38.all.events.ioe
suppa.py参数说明：（通过conda安装的话，不需要再前面加python来运行命令的）
generateEvents suppa.py子命令用于生成可变剪切时间，后面的参数是该子命令参数
-i：gtf注释文件
-o 输出文件名前缀
-f 可变剪切注释格式，分: ioe for local events, ioi for transcript events.
-e : 只用于ioe格式， 有这些事件
SE: Skipping exon (SE) events
SS: Alternative 5' (A5) and 3' (A3) splice sites (it generates both)
MX: Mutually Exclusive (MX) exons
RI: Retained intron (RI)
FL: Alternative first (AF) and last (AL) exons (it generates both)
计算local AS event 的psi value, 需要转录本tpm
获取转录本tpm
mkdir ../as
cd ../as
multipleFieldSelection.py -i ~/Salmon_output/*/quant.sf -k 1 -f 4 -o iso_tpm.txt
head iso_tpm.txt
HMDM_M0_rep1_SRR11652531_output HMDM_M0_rep2_SRR11652533_output HMDM_M0_rep3_SRR11652535_output HMDM_M1_rep1_SRR11652532_output HMDM_M1_rep2_SRR11652534_output HMDM_M1_rep3_SRR11652536_output
ENST00000456328.2 0.000000 0.047879 0.223474 0.000000 0.000000 0.000000
ENST00000450305.2 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000488147.1 3.598279 5.462875 3.451416 3.008844 6.047462 4.102420
ENST00000619216.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000473358.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000469289.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000607096.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000417324.1 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
ENST00000461467.1 0.000000 0.000000 0.000000 0.088777 0.000000 0.000000
multipleFieldSelection.py参数说明：
-i ：输入文件，使用了通配符*来匹配多个文件
-k : 指输入文件的共用的字段，这里设为1，指quant.sf的第一列的转录本ID是所有文件共用的id。
-f: 提取第4 列值，第四列是tpm值。
-o：输出文件名
计算AS events 的psi值。
suppa.py psiPerEvent -i ../ref/gencode.v38.all.events.ioe -e iso_tpm.txt -o HMDM_as_event 1>psiPerEvent_log.txt 2>&1 &
head HMDM_as_event.psi
HMDM_M0_rep1_SRR11652531_output HMDM_M0_rep2_SRR11652533_output HMDM_M0_rep3_SRR11652535_output HMDM_M1_rep1_SRR11652532_output HMDM_M1_rep2_SRR11652534_output HMDM_M1_rep3_SRR11652536_output
ENSG00000000003.15;A5:chrX:100635746-100636191:100635746-100636608:- 0.0 0.0 0.0 0.0 0.0 0.0
ENSG00000000003.15;A5:chrX:100635746-100636608:100635746-100636793:- 0.6043555603960755 1.0 1.0 1.0 1.0 1.0
ENSG00000000003.15;AF:chrX:100635746-100636191:100636689:100635746-100636793:100637104:- 0.0 nan nan nan nan nan
ENSG00000000003.15;SE:chrX:100630866-100632485:100632568-100633405:- 1.0 1.0 1.0 1.0 1.0 1.0
ENSG00000000419.14;A3:chr20:50940955-50941105:50940933-50941105:- 0.0 0.0 0.03959799699675568 0.0 0.1305081076991619 0.0
ENSG00000000419.14;A3:chr20:50940955-50942031:50940933-50942031:- 0.02079990366574006 0.009161120464883297 0.0051864495628405225 0.017888399727889845 0.014109994702014904 0.03745825846097759
参数说明：
-i : ioe 文件
-e : 转录本表达量
-o 输出文件名
输出的内容的第一列，是local AS event ID，剩下的列是psi value.
运行过程可能会产生类似下面这种ERROR的错误，这是由于没有在tpm表达量矩阵中找到该转录本。
ERROR:psiCalculator:PSI not calculated for event ENSG00000177873.14;SE:chr3:40482704-40483653:40483723-40486806:+.
ERROR:psiCalculator:transcript ENST00000640853.1 not found in the "expression file".
ERROR:psiCalculator:PSI not calculated for event ENSG00000145012.14;SE:chr3:188225527-188312835:188312918-188406112:+.
ERROR:psiCalculator:transcript ENST00000640853.1 not found in the "expression file".
我们在psiPerEvent_log.txt 日志里找下有多少转录本没有找到
grep "not found" psiPerEvent_log.txt | wc -l
1264
有个1264转录本没找到，我们在看下基因注释文件和iso_tpm文件的转录本有多少。
wc -l ../salmon_output/*/quant.sf
236187 ../salmon_output/HMDM_M0_rep1_SRR11652531_output/quant.sf
236187 ../salmon_output/HMDM_M0_rep2_SRR11652533_output/quant.sf
236187 ../salmon_output/HMDM_M0_rep3_SRR11652535_output/quant.sf
236187 ../salmon_output/HMDM_M1_rep1_SRR11652532_output/quant.sf
236187 ../salmon_output/HMDM_M1_rep2_SRR11652534_output/quant.sf
236187 ../salmon_output/HMDM_M1_rep3_SRR11652536_output/quant.sf
1417122 total
wc iso_tpm.txt -l
236187 iso_tpm.txt
cut -f 3 ../ref/gencode.v38.annotation.gtf | grep "transcript" | wc -l
237012
237012 - 236187 = 825 ，可以知道salmon定量出来，也是有一些转录本没定量出来的。个人觉得，suppa2报错说找不到的转录本，并不多，应该问题不大。虽然还有部分转录本不知去哪了。
注：若是根据我的笔记，从gencode 上下载的序列文件，以及按照代码来，iso_tpm.txt 文件是可以直接使用的。iso_tpm.txt 里的ID已经只是ENSEMBLE的转录本ID了。这是因为我在构建salmon index时加了--gencode，程序自动处理了ID。若是那步没加，或者使用了其他的参考序列。可根据实际进行处理
参考生信技能树代码：
perl -alne '{/(\|.*\|)\t/; ;s/$1//g;s/\|//g;print}' iso_tpm.txt > iso_tpm_formatted.txt
做个箱型图，查看某个AS event的psi value.
## 官方提供的脚本，
wget https://raw.githubusercontent.com/comprna/SUPPA/master/scripts/generate_boxplot_event.py
python generate_boxplot_event.py -i HMDM_as_event.psi -e "ENSG00000004866.22;AF:chr7:116953327:116953691-117099762:117014888:117015014-117099762:+" -g 1-3,4-6 -c M0,M1 -o .
参数说明：
-i 输入psi文件
-e AS event id
-g 分组， 1-3 为一组，4-6为一组
-c 分组命名
-o 输出目录
差异可变剪切事件寻找, 这里需要对tpm,和psi进行分组。
### 因为它的第一列，是id,但是没有id 列名，被样本名给占了。
### 虽然按照之前生信技能树文章代码里的注释说过表头似乎是并不是很重要
### 这里还是保留表头吧
awk -v OFS='\t' '{if (NR == 1) {print $1,$2,$3} else {print $1,$2,$3,$4}}' HMDM_as_event.psi > M0.psi
awk -v OFS='\t' '{if (NR == 1) {print $4,$5,$6} else {print $1,$5,$6,$7}}' HMDM_as_event.psi > M1.psi
awk -v OFS='\t' '{if (NR == 1) {print $1,$2,$3} else {print $1,$2,$3,$4}}' iso_tpm.txt > M0.tpm
awk -v OFS='\t' '{if (NR == 1) {print $4,$5,$6} else {print $1,$5,$6,$7}}' iso_tpm.txt > M1.tpm
suppa.py diffSplice -m empirical -gc -i ../ref/gencode.v38.all.events.ioe -p M0.psi M1.psi -e M0.tpm M1.tpm -o HMDM_ds
ll -h HMDM_ds* |cut -d" " -f 5-
24M Jun 26 17:15 HMDM_ds.dpsi
31M Jun 26 17:15 HMDM_ds.psivec
head HMDM_ds.dpsi
M00-M11_dPSI M00-M11_p-val
ENSG00000000003.15;A5:chrX:100635746-100636191:100635746-100636608:- 0.0 0.1513486513
ENSG00000000003.15;A5:chrX:100635746-100636608:100635746-100636793:- 0.1318814799 0.1513486513
ENSG00000000003.15;AF:chrX:100635746-100636191:100636689:100635746-100636793:100637104:- nan 1.0
ENSG00000000003.15;SE:chrX:100630866-100632485:100632568-100633405:- 0.0 0.1513486513
ENSG00000000419.14;A3:chr20:50940955-50941105:50940933-50941105:- 0.0303033702 0.4025065843
ENSG00000000419.14;A3:chr20:50940955-50942031:50940933-50942031:- 0.0114363931 0.4025065843
ENSG00000000419.14;A5:chr20:50940933-50941105:50940933-50941129:- 0.0461316111 0.4025065843
ENSG00000000419.14;A5:chr20:50942126-50945737:50942126-50945847:- -0.0835089009 0.4025065843
ENSG00000000419.14;AF:chr20:50955285-50958363:50958417:50955285-50958741:50959140:- nan 1.0
diffSplice 参数：
-m calculate the significance (empirical/classical)
-gc Correction of the p-values by gene.
-i Input file with the local or transcript events, .ioe or .ioi format, respectively
-p PSI files.
-e Transcript expression (TPM) files.
-o Name of the output
查看下具有统计学意义上的结果
cat HMDM_ds.dpsi|perl -alne '{print if $F[2] <0.01}' |wc -l
55
cat HMDM_ds.dpsi|perl -alne '{print if $F[2] <0.05}' |wc -l
532
cat HMDM_ds.dpsi|perl -alne '{print if $F[2] <0.05}' | head
M0-M1_dPSI M0-M1_p-val
ENSG00000001167.15;SE:chr6:41079164-41080811:41080897-41084046:+ -0.3118838972 0.0154845155
ENSG00000002587.10;AF:chr4:11400113-11428463:11428597:11400113-11428699:11428894:- 0.3332659782 0.042957043
ENSG00000002587.10;AF:chr4:11400113-11428463:11428597:11400113-11429425:11429564:- 0.6431379821 0.042957043
ENSG00000002587.10;AF:chr4:11400113-11428699:11428894:11400113-11429425:11429564:- 0.3286477452 0.042957043
ENSG00000003400.15;SE:chr2:201195948-201203730:201203766-201208075:+ -0.4890349374 0.044955045
ENSG00000004468.13;A5:chr4:15825016-15834217:15824998-15834217:+ -0.3709872401 0.0157342657
ENSG00000004468.13;SE:chr4:15816640-15824881:15824998-15834217:+ 0.100261049 0.0404595405
ENSG00000004468.13;SE:chr4:15816640-15824881:15825016-15834217:+ -0.3018835422 0.0157342657
ENSG00000006451.8;SE:chr7:39696859-39697408:39697510-39706123:+ 0.3106216186 0.012987013DTU，是同一基因内，在不同条件下，转录本表达量占比发生变化，如下图3,在condition1 和condition2条件下，GeneA,的转录本1，2的表达量占比发生变化，是DTU。转录本3不是。所以如果一个基因只有一个转录本时，表达量占比不会发生改变，那它就不是DTU了。
DTU分析和之前local AS event一样还是分为三步。
生成transcript event
suppa.py generateEvents -f ioi -i ../ref/gencode.v38.annotation.gtf -o ../ref/gencode.v38.isoforms
head ../ref/gencode.v38.isoforms.ioi
seqname gene_id isoform_id inclusion_transcripts total_transcripts
chr1 ENSG00000223972.5 ENSG00000223972.5;ENST00000456328.2 ENST00000456328.2 ENST00000456328.2,ENST00000450305.2
chr1 ENSG00000223972.5 ENSG00000223972.5;ENST00000450305.2 ENST00000450305.2 ENST00000456328.2,ENST00000450305.2
chr1 ENSG00000243485.5 ENSG00000243485.5;ENST00000473358.1 ENST00000473358.1 ENST00000473358.1,ENST00000469289.1
chr1 ENSG00000243485.5 ENSG00000243485.5;ENST00000469289.1 ENST00000469289.1 ENST00000473358.1,ENST00000469289.1
chr1 ENSG00000284332.1 ENSG00000284332.1;ENST00000607096.1 ENST00000607096.1 ENST00000607096.1
chr1 ENSG00000268020.3 ENSG00000268020.3;ENST00000606857.1 ENST00000606857.1 ENST00000606857.1
chr1 ENSG00000240361.2 ENSG00000240361.2;ENST00000642116.1 ENST00000642116.1 ENST00000642116.1,ENST00000492842.2
chr1 ENSG00000240361.2 ENSG00000240361.2;ENST00000492842.2 ENST00000492842.2 ENST00000642116.1,ENST00000492842.2
chr1 ENSG00000186092.7 ENSG00000186092.7;ENST00000641515.2 ENST00000641515.2 ENST00000641515.2
计算transcript 的psi, 这里的psi其实就是某个转录本tpm值占所有转录本tpm值和的比值。
suppa.py psiPerIsoform -g ../ref/gencode.v38.annotation.gtf -e iso_tpm.txt -o iso
head iso_isoform.psi
HMDM_M0_rep1_SRR11652531_output HMDM_M0_rep2_SRR11652533_output HMDM_M0_rep3_SRR11652535_output HMDM_M1_rep1_SRR11652532_output HMDM_M1_rep2_SRR11652534_output HMDM_M1_rep3_SRR11652536_output
ENSG00000223972.5;ENST00000456328.2 nan 1.0 1.0 nan nan nan
ENSG00000223972.5;ENST00000450305.2 nan 0.0 0.0 nan nan nan
ENSG00000243485.5;ENST00000473358.1 nan nan nan nan nan nan
ENSG00000243485.5;ENST00000469289.1 nan nan nan nan nan nan
ENSG00000284332.1;ENST00000607096.1 nan nan nan nan nan nan
ENSG00000268020.3;ENST00000606857.1 nan nan nan nan nan nan
ENSG00000240361.2;ENST00000642116.1 nan nan nan nan nan nan
ENSG00000240361.2;ENST00000492842.2 nan nan nan nan nan nan
ENSG00000186092.7;ENST00000641515.2 nan nan nan nan nan nan
nan是由于转录本tpm值为0. 当了分母了。
DTU计算, 这里需要对tpm,和psi进行分组，tpm之前分好了，这里分下iso_isoform.psi。
awk -v OFS='\t' '{if (NR == 1) {print $1,$2,$3} else {print $1,$2,$3,$4}}' iso_isoform.psi > M0_iso.psi
awk -v OFS='\t' '{if (NR == 1) {print $4,$5,$6} else {print $1,$5,$6,$7}}' iso_isoform.psi > M1_iso.psi
suppa.py diffSplice -m empirical -gc -i ../ref/gencode.v38.isoforms.ioi -p M0_iso.psi M1_iso.psi -e M0.tpm M1.tpm -o HMDM_ds_iso
cat HMDM_ds_iso.dpsi|perl -alne '{print if $F[2] <0.01}' |wc -l
134
cat HMDM_ds_iso.dpsi|perl -alne '{print if $F[2] <0.05}' |wc -l
1302
cat HMDM_ds_iso.dpsi|perl -alne '{print if $F[2] <0.05}' | head
M0_iso-M1_iso_dPSI M0_iso-M1_iso_p-val
ENSG00000000457.14;ENST00000367772.8 0.2801505809 0.04995005
ENSG00000001167.15;ENST00000341376.11 -0.3118838972 0.0074925075
ENSG00000001167.15;ENST00000353205.5 0.3118838972 0.0074925075
ENSG00000002587.10;ENST00000510712.1 -0.3222232746 0.014985015
ENSG00000002587.10;ENST00000514690.5 0.3454912316 0.014985015
ENSG00000004468.13;ENST00000226279.8 -0.2746099506 0.0237262737
ENSG00000004468.13;ENST00000502843.5 0.2018218955 0.0237262737
ENSG00000004468.13;ENST00000510674.1 0.0797545927 0.0457875458
ENSG00000004660.15;ENST00000348335.7 -0.4477146765 0.035964036
转录组练习之测试流程和代码，jimmy
https://github.com/COMBINE-lab/salmon
https://combine-lab.github.io/alevin-tutorial/2019/selective-alignment/
https://github.com/comprna/SUPPA/wiki/SUPPA2-tutorial
http://yulab-smu.top/clusterProfiler-book/chapter11.html
https://mp.weixin.qq.com/s/kG2uI3me8Xo69mClThRutA
https://mp.weixin.qq.com/s/NC5mKwAsYbm1Q-PFC5IEOg
https://mp.weixin.qq.com/s/YjraGRBCc725IfhYB1VPKQ