Bulk RNA-seq上游分析流程（Salmon）

转录组分析一般通过比对或非比对的方法对基因表达进行定量，得到一个样本-基因表达矩阵用于下游数据分析。下游分析通常是进行差异表达分析，对差异基因进行功能注视，如Gene Ontology (GO) 或 KEGG pathway。也可以对差异基因进行蛋白互作网络 (PPI) 分析，找出可能存在的关键基因。除了差异表达分析，还可以进行WGCNA分析，多与临床数据相关。

原始测序数据质控

测序数据有可能含有测序接头和质量较低的reads，去除这些序列对于后续分析非常重要，因此先用fastp软件进行质控，去除测序接头和低质量reads。fastp的优势在于不用输入接头类型，它会自动进行检测，接头去除的比较干净，另外对于双端测序还可以进行碱基校正。

fastp -l 25 -w 4 --detect_adapter_for_pe -i input_R1.fastq -I input_R2.fastq -o output_R1.fastq -O output_R2.fastq

参数解释：

-l 表示保留最终序列长度大于25的reads
-w 表示线程
--detect_adapter_for_pe 双端测序参数
-i, -I 分别表示read1，read2两个输入文件
-o, -O 分别表示read1，read2对应的输出文件

利用salmon进行基因表达定量

salmon是一款不通过序列比对就可以快速完成生物学定量的RNA-seq数据分析工具。它的使用流程包括两步：1.建立索引 2.对reads进行基因表达定量(quantification)。

第一步，提取转录本-基因对应关系

转录谱文件可以到Ensembl下载, 转录谱参考文件内容如下：

>ENST00000632684.1 cdna chromosome:GRCh38:7:142786213:142786224:1 gene:ENSG00000282431.1 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene gene_symbol:TRBD1 description:T cell receptor beta diversity 1 [Source:HGNC Symbol;Acc:HGNC:12158]
GGGACAGGGGGC

第一行是header信息，以>开头，第一列 ENST00000632684.1 为转录本ID，salmon默认对转录本进行定量，但是一个基因有多个转录本，因此需要把同一基因的各个转录本定量结果合并在一起，所以需要提取一个转录本-基因对应的信息文件。gene_symbol:TRBD1中的TRBD1即为基因名字，把第一列和基因列提取出来就行了。

提取命令如下：

zgrep ">" Homo_sapiens.GRCh38.cdna.all.fa.gz | sed 's/>//g' | sed 's/cdna.*gene_symbol://g' | sed 's/description.*//g' > gene_map.txt

注意：有些物种的EnsembI转录组参考序列只有转录本ID（ENST00000632684.1）和基因ID（ENSG00000282431.1），部分基因名字是缺失的，提取时需留意！！！无基因名字的可以用基因ID代替

第二步，建立索引

salmon index -p 4 -t Homo_sapiens.GRCh38.cdna.all.fa.gz -i salmon_index -k 31 --type quasi

参数解释：

-p 表示线程
-t 转录谱参考序列，可以到ensembl下载
-i 索引文件输出的文件夹
-k 索引k-mers长度
--type 索引类型

第三步，基因定量

salmon quant -p 4 --validateMappings -i salmon_index -l A -g gene_map.txt -1 output_R1.fastq -2 output_R2.fastq -o Salmon_out/output_name

参数解释：

-p 表示线程
--validateMappings 启用选择性比对到转录本。能够提高比对的灵敏度和特异性从而提高了定量的准确度。
-l library 类型，默认A，自动检测
-g 转录本-基因之间的对应关系，如果不指定，不会输出对基因表达的定量，而是输出对转录本表达的定量
-1,-2 质控后的read1，read2两个fastq文件
-o 输出文件夹，Salmon_out是所有样本输出的目录，output_name是本次running的样本定量结果输出目录

输出目录结构如下：

Salmon_out
├── CC_Cas9-1
├── CC_Cas9-2
├── CC_Cas9-3
├── CC_Cas9_Cre-1
├── CC_Cas9_Cre-2
├── CC_Cas9_Cre-3
├── CC_Cas9_Flpo-1
├── CC_Cas9_Flpo-2
└── CC_Cas9_Flpo-3

每个目录下文件结构如下：

Salmon_out/CC_Cas9-1
├── aux_info
├── cmd_info.json
├── libParams
├── lib_format_counts.json
├── logs
├── quant.genes.sf
└── quant.sf

其中quant.sf为转录本定量结果：

Name    Length  EffectiveLength TPM     NumReads
ENST00000474626.5       639     522.973 75.863174       32.032
ENST00000476222.5       685     690.231 8.012039        4.465
ENST00000496813.5       642     654.652 2.670470        1.411
ENST00000490191.5       996     844.000 0.000000        0.000
ENST00000472218.1       1109    957.000 0.000000        0.000
ENST00000479802.1       413     261.136 281.314078      59.311

如果指定了参数-g，则会输出基因定量结果，存储在文件quant.genes.sf中：

Name    Length  EffectiveLength TPM     NumReads
AC008626.1      255     775.152 202.283 126.596
VN1R107P        209     80      0       0
Z97832.1        328     179     0       0
OR7E13P 910     554.755 2.23267 1
AC006511.1      315     167     0       0
MTCO3P12        547     355.804 8.42352 2.42
MTCO2P12        682     531     0       0
MTND2P28        1044    730.596 29.969  17.678

利用自定义脚本生成表达矩阵

用shell脚本提取基因名称和基因定量两列数据，以table分割

#!/bin/bash
# 脚本名称：extract_gene_cols.sh

for i in `ls $1`
do
    sample_name=$i
    awk '{print$1,$5}' $1/$i//quant.genes.sf | grep -v 'Name' | sed 's/ /\t/' > $1/$i/${i}.count
done

使用方法: bash extract_gene_cols.sh Salmon_out, 参数Salmon_out即为所有样本的输出文件夹。

用python脚本合并创建表达矩阵

#!/usr/local/bin/python
# -*- coding:utf-8 -*-
# 脚本名：generate_count_matrix.py
import sys
import pandas as pd
from glob import iglob

sample_path = iglob(sys.argv[1] + "/*/*.count")
dfs = []

for df in sample_path:
    fn = df.split('/')[-1]
    sample_name = fn[:-6]
    dfs.append(pd.read_table(df, header=None, index_col=0,sep='\t', names=["geneid", sample_name]))

matrix = pd.concat(dfs, axis=1)
matrix.to_csv(sys.argv[2])

使用方法: python generate_count_matrix.py Salmon_out Count_matrix.csv, 第一个参数Salmon_out即为所有样本的输出文件夹, 第二个参数Count_matrix.csv为输出表达矩阵文件名字。

最终生成的Count_matrix.csv表达矩阵为行为基因，列为样本的表达矩阵：

geneid,CC_Cas9_Flpo-1,CC_Cas9-1,CC_Cas9_Flpo-2,CC_Cas9_Cre-1,CC_Cas9-2,CC_Cas9_Cre-3,CC_Cas9_Flpo-3,CC_Cas9-3,CC_Cas9_Cre-2
AC217785.3,2.639,0.0,1.9909999999999999,0.0,0.0,0.0,0.0,0.0,0.0
FRG1KP,0.0,0.0,0.0,0.0,0.0,1.0,0.0,1.915,0.0
AC007322.4,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0
RPS6P20,511.55699999999996,975.587,1573.53,1309.35,689.59,1340.45,308.182,752.914,217.21400000000003
OR13I1P,0.986,6.224,20.303,5.135,4.225,10.315,0.0,3.5069999999999997,2.539