ESPRESSO全称：Error Statistics PRomoted Evaluator of Splice Site Options，一种处理long read RNA-seq数据alignment的方法，可以有效地提高splice junction精度和isoform定量。
下载地址：https://github.com/Xinglab/espresso/releases/
省流版
用基因组+参考注释确定高置信度的splice junction，根据这些SJ对基因组中的每个reads进行校正和恢复，最后用校正过的SJ和注释信息对所有isoforms定量，输出更新后的注释信息gtf文件和isoforms表达丰度文件。
详细版
依赖的软件：
- Perl >= 5.8 built with threading enabled
- hmmer >= 3.3.1
conda install -c bioconda hmmer
- BLAST >= 2.8.1
软件包地址：https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/ - samtools >= 1.6
（最好创建一个新的环境来安装samtools）
conda create -n samtools_env
conda activate samtools_env
conda install -c bioconda samtools openssl=1.0
（以下是可视化依赖的软件）
- UCSC tools
-- bedGraphToBigWig
-- faToTwoBit
-- twoBitInfo
软件包地址：https://hgdownload.soe.ucsc.edu/admin/exe/
可以用rsync -aP rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/ ./下载，默认是当前目录，但是有n个执行文件，最好创建一个新的目录统一存放，不然当前目录会爆炸 - Python 3
-- NumPy - IGV：https://igv.org/
需要的文件：
- sam文件/bam文件
如果不是这种文件类型，可以用minimap2或ONT guppy basecaller输出这种文件 - gtf注释文件
- fasta参考基因组文件
使ESPRESSO在所有目录下可用
- 在每个pl或py脚本文件第一行加
#! /usr/bin/perl或#! /usr/bin/python3 - 更改每一个脚本文件权限，增加可执行权限
chmod +x ESPRESSO_C.pl - 将ESPRESSO执行文件加入环境变量
这种方法可以在本用户下永久添加环境变量
vim ~/.bashrc
输入
export PATH=/data/workdir/wubq/soft/espresso_v_1_3_2/src/:$PATH
export PATH=/data/workdir/wubq/soft/espresso_v_1_3_2/visualization/:$PATH
保存
source ~/.bashrc
如果你没有做这步的话，后续用pl或py脚本需要在命令行前加perl或python3，并且执行文件要加路径
预备工作：准备samples.tsv文件
/path/to/SIRV2_3.sort.sam test_sample_name # tab键分隔
在本例中，只有一个SAM文件。如果有多个SAM文件，那么每个文件都单独一行。第二列是sample名称，可用于在最终输出中将多个输入SAM文件分组在一起。
第一步：确定高置信度的splice junctions——ESPRESSO_S
用法：
ESPRESSO_S.pl -L samples.tsv -F ref.fa -A anno.gtf -O work_dir
-L 刚刚创建的tsv文件，第一列是带有绝对路径的sorted sam文件；第二列是sample名称
-F reference文件
-A 注释文件
-B 你认为信得过的splice junction，bed文件格式
-O 输出文件的目录名称（不是输出文件名），默认在当前目录
-T 线程数
用test_sirv数据跑的话输出是：
[Thu Sep 7 10:42:45 2023] Loading reference
Worker 0 begins to scan:
/data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_data_espresso_sirv/SIRV2_3.sort.sam
Worker 0 finished reporting.
[Thu Sep 7 10:42:50 2023] Re-cluster all reads
[Thu Sep 7 10:42:50 2023] Loading annotation
[Thu Sep 7 10:42:50 2023] Summarizing annotated splice junctions for each read group
/data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_data_espresso_sirv/SIRV2_3.sort.sam(0)
0_1(0)
Worker 0 begins to scan:
/data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_data_espresso_sirv/SIRV2_3.sort.sam
/data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_data_espresso_sirv/SIRV2_3.sort.sam 0
Worker 0 finished reporting.
[Thu Sep 7 10:42:57 2023] ESPRESSO_S finished its work.
在test_sirv目录下有以下文件：
ESPRESSO_S结果
其中sample.tsv.updated文件是为SAM文件分配目标id
里面只分配到一个，ID为0
/path/to/SIRV2_3.sort.sam test_sample_name 0
第二步：校正和恢复每个read的splice junction——ESPRESSO_C
用法：
ESPRESSO_C.pl -I test_sirv -F SIRV2.fasta -X 0 -T 5
-I (in)由上一步ESPRESSO_S所得的文件所在目录
-F fasta reference文件
-X 目标ID，如果有多个输入，那么ESPRESSO_C将为samples.tsv.updated中的每个ID分别运行
-T 线程数
输出：
[Thu Sep 7 13:47:51 2023] Loading splice junction info
[Thu Sep 7 13:47:52 2023] Requesting system to split SAMLIST into 5 pieces
Divided SAM(LIST) sizes:
sam.list3aa 1242019
sam.list3ab 1242019
sam.list3ac 1242019
sam.list3ad 1242019
sam.list3ae 1242019
SAM(LIST) was divided successfully.
First group of divided SAM(LIST) files:
sam.list3aa: 0
First reads were recorded successfully for 1 files.
[Thu Sep 7 13:47:52 2023] Loading references
[Thu Sep 7 13:47:52 2023] Scanning SAMLIST by 2 workers
Worker 1 begins to scan sam.list3aa.
0 16
# 此处省略一堆
Building a new DB, current time: 09/07/2023 13:48:14
New DB name: /data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_sirv/0/blast_0/current_db
New DB title: current_db
Sequence type: Nucleotide
Deleted existing Nucleotide BLAST database named /data/workdir/wubq/soft/espresso_v_1_3_2/test_data/test_sirv/0/blast_0/current_db
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 2 sequences in 0.000781059 seconds.
Worker 1 finished reporting.
[Thu Sep 7 13:48:37 2023] ESPRESSO_C finished its work.
第三步：从reads中鉴定和定量所有isoforms——ESPRESSO_Q
用法：
ESPRESSO_Q.pl -A SIRV_C.gtf -L test_sirv/samples.tsv.updated -V test_sirv/samples_N2_R0_compatible_isoform.tsv
-A gtf注释文件
-L sample list的tsv文件
-V 对于每个reads的兼容isoform的tsv输出文件（可能是从属于某个isoform的read的可变剪切类型）
-O 输出文件所在目录，默认在sample list的tsv文件目录里
-T 线程数
输出：
[Thu Sep 7 14:02:47 2023] Loading annotation
[Thu Sep 7 14:02:47 2023] Summarizing annotated isoforms
[Thu Sep 7 14:02:47 2023] Loading corrected splice junctions and alignment information by ESPRESSO
[Thu Sep 7 14:02:48 2023] Categorizing reads according to annotation
[Thu Sep 7 14:02:48 2023] thread_id: 4 starting to process chr: SIRV2
running command: /usr/bin/perl ../src/ESPRESSO_Q_Thread.pm test_sirv/thread_shared_arguments.tmp test_sirv/thread_SIRV2_arguments.tmp
finished command: /usr/bin/perl ../src/ESPRESSO_Q_Thread.pm test_sirv/thread_shared_arguments.tmp test_sirv/thread_SIRV2_arguments.tmp
[Thu Sep 7 14:02:58 2023] thread_id: 4 finished
[Thu Sep 7 14:02:58 2023] cleaning up threads
[Thu Sep 7 14:03:00 2023] ESPRESSO finished quantification
第四步：可视化
用法：
visualize.py --genome-fasta SIRV2.fasta --updated-gtf test_sirv/samples_N2_R0_updated.gtf --abundance-esp test_sirv/samples_N2_R0_abundance.esp --target-gene SIRV2 --minimum-count 1 --descriptive-name SIRV --output-dir test_sirv/visualization
--genome-fasta 输入参考基因组文件
--updated-gtf 由ESPRESSO生成的，更新后的gtf文件
--abundance-esp 由ESPRESSO生成的丰度文件
--target-gene 要可视化的基因，使用gene_id来匹配ESPRESSO输出的新异构体
--minimum-count 只考虑在一个样本的count满足最小count的isoforms (我不理解这是什么意思TAT)
(原句：only isoforms where the (normalized) count for a sample meets the minimum count are considered)
--normalize-counts-to-cpm 将原始计数转换为CPM
--descriptive-name 图标及文件名
--output-dir 输出文件所在目录
可视化输出文件
这些输出文件需要用IGV可视化
- Genomes -> Load Genome from File -> {file}.fasta and {file}.fasta.fai
- File -> Load from File -> {file}.gtf
-- Right click gtf track.
-- Select "squished".
-- Set track height to a large enough value. - File -> Load from File -> {sample_name}.bw for each sample
-- Right click track.
-- Select "Autoscale".
--- Select "Maximum" under "Windowing Function".
-- "Change Track Color". - View -> Add New Panel
-- Click and drag panel borders to show the new panel. - File -> Load from File -> target_genes/{file}.bed
-- Click name of new track and drag to the appropriate panel. - Adjust the coordinates to show the data.
- File -> Save Image ->
-- Change file extension to .svg
visualization_sirv.png
解释一下为什么转换的是CPM，而不是平常转录组分析中表达定量常用的RPKM、FPKM、TPM
CPM全称是Counts per million，是对测序深度进行标准化，而TPM是先对基因长度进行标准化，再对测序深度进行标准化。short read测序建库时切割了很多个比较均匀的短read，基因表达越多，reads数越多，所以需要校正。而long read测序，一个read几乎就是一条transcript，所以不需要校正基因长度。3‘端polyA如果出错，5’端可能会测不到或不均匀，也没有校正的意义。ESPRESSO用的数据是三代测序技术的数据。
推荐另一个可视化软件，也是邢毅老师团队开发的，rMATS-long：https://github.com/Xinglab/rMATS-long。可视化效果更好。但是它有一点不太好的地方是只能画丰度最高的5个或以下isoform，其他一律归为”others“。可以在visualize_isoforms.py中增加颜色，增加越多越好，增加多一种颜色就可以画多一条isoform。在画图时通过
max-transcripts参数输入你的isoform个数，即可全部画出。
rMATS-long可视化isoform占比和CPM
rMATS-long可视化isoform结构
visualize_isoforms.py增加颜色
解释一下各输出文件
sample_N2_R0_updated.gtf：为检测到的isoforms更新的GTF注释文件
- 每一个检测到的isoform会视为一个转录本
- source列 表示每个异构体是“novel_isoform”还是“annotated_isoform”
sample_N2_R0_abundance.esp：检测到的isoforms的表达丰度文件，是tsv文件
- 每个检测到的isoform会在单独的行上写明
- 第一列是 transcript_ID, transcript_name, gene_ID
- 在samples.tsv中提供的每个样本名称都有一个额外的列。这些列表示从该样本中被算到这个isoform的reads数。
- counts是通过期望最大化（EM）来分配的。每个输入read最多贡献1个count，可以是单个isoform，也可以是多个isoform。
sample_N2_R0_compatible_isoform.tsv：每个read的兼容isoform（可能是从属于某个isoform的read？猜的）
- 该文件仅在使用ESPRESSO_Q的-V参数时生成。
- 各列分别是 read_id, sample_name, read_classification, compatible_isoforms
- 有以下几种分类 NIC/NNC, NCD, ISM, FSM, single-exon
####################################
2023.12.1更新
报错合集：
- ESPRESSO_S报错：
- out of memory
原因：1.3.2版本及以上需要perl storable版本在3.0以上，我们的服务器是2.62。可以通过perl -e 'use Storable; print("gstorable::VERSIONn")!查看版本
解决：升级perl storable
out of memory
- ESPRESS_C报错：
- 找不到3.3.1版本的nhmmer，但实际上安装了nhmmer，版本号更高
原因：1.3.2版本设定只接受3.3.1的nhmmer
解决：自己改一下ESPRESSO_C.pl脚本，把==1改为>=1，如此类推。或者直接用最新版本1.4
nhmmer
源代码
- ESPRESSO_Q报错：
- 没有一个git库(not a git repository)
原因：未知
解决：用git init初始化后也没用，还是不行。
不存在git库
