转录组上游分析流程(四）

环境部署——数据下载——查看数据(非质控)——数据质控——数据过滤(过滤低质量数据)——数据比对及定量

数据比对：

1、参考基因组准备：

Ensembl官网

左上箭头分别是最新版本号和Fasta文件下载链接

下载primary_assembly.fa.gz文件(复制网页链接+下边需要下载的内容)

2.基因组注释文件准备：Ensembl官网

右侧有Gene annotation模块。

Download FASTA 可以下载fasta文件，包含了gene和cDNA(转录组)；

Download GTF or GFF3 files for genes，cDNAs，ncRNA，proteins包含了参考基因组注释文件，gtf个gff 下载注释文件的时候需要注意版本信息，要前后一致。

比如wget -c ftp://ftp.ensembl.org/pub/release 113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

关键就是所有的release信息需要对应起来

## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# ftp://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/

# 下载基因组序列
wget -c ftp://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# 下载基因组注释文件
wget -c ftp://ftp.ensembl.org/pub/release-113/gtf/homo_sapiens/Homo_sapiens.GRCh38.113.gtf.gz
wget -c ftp://ftp.ensembl.org/pub/release-113/gtf/homo_sapiens/Homo_sapiens.GRCh38.113.chr.gtf.gz

# 后台运行 注意版本号
nohup wget -c ftp://ftp.ensembl.org/pub/release-113/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz > dan.log &
nohup wget -c ftp://ftp.ensembl.org/pub/release-113/gtf/homo_sapiens/Homo_sapiens.GRCh38.113.gtf.gz > dan.log &

# 检查大文件完整性
gzip -t *.gz

3.fastq与fasta文件转换：

转换成fasta的目的是去除附加和质量控制信息，便于后续分析。

方法1：

1. zless -S ./trim_galore/SRR23881762_1_val_1.fq: zless：用于查看压缩文件内容的命令，类似 less，但它支持 .gz 压缩文件。-S：使 zless 水平滚动，不会自动换行。./trim_galore/SRR23881762_1_val_1.fq：文件路径，表示查看 SRR23881762_1_val_1.fq 文件。
2. | awk '{ if(NR%4==1){print">" substr(0,2)} if(NR%4==2){print} }': awk：文本处理工具，按行处理数据。0：在 awk 中表示当前行的整个内容。NR%4==1：表示每4行中第1行，因为 FASTQ 文件中每个序列都是4行组成的（@序列ID、序列、+、质量分值），所以第1行是序列ID行。print ">" substr($0,2)：将 @ 开头的序列ID行替换成 > 开头，并从第二个字符开始显示（即去掉原来的 @）。NR%4==2：表示每4行中的第2行，这一行是实际的序列内容，所以直接打印。
3. | less -S: less：分页查看工具。-S：同样启用水平滚动。

方法2：

1. zless -S ./trim_galore/SRR23881762_1_val_1.fq: 与上面的解释相同，用 zless 查看压缩的 FASTQ 文件内容，并启用水平滚动。
2. | paste - - - -: paste：用于将多行合并成一行的命令。-：每次读取4行，合并成一行（用 TAB 分隔）。这一步的作用是将 FASTQ 文件中的每个序列（4行）合并成一行。
3. | cut -f 1,2: cut：用于从文本中提取指定字段的命令。-f 1,2：表示提取合并后的第1和第2个字段，第1字段是序列ID（原来的第1行），第2字段是序列内容（原来的第2行）。
4. | tr '@' '>': tr：用于替换或删除字符的命令。'@' '>'：将序列ID中的 @ 替换为 >，符合 FASTA 格式的要求。
5. | tr '\t' '\n': tr '\t' '\n'：将 TAB 替换为换行，将原来 paste 合并的一行再次拆分为两行（序列ID和序列）。
6. | less -S: 分页查看最终结果。

# 方法1
zless -S ./trim_galore/SRR23881762_1_val_1.fq |awk '{ if(NR%4==1){print">" substr($0,2)} if(NR%4==2){print}  }' | less -S
# Tips：substr($1, 2)，从第一个字段里的第2个字符开始，一直到结束

# 方法2
zless -S ./trim_galore/SRR23881762_1_val_1.fq |paste - - - - |cut -f 1,2 |tr '@' '>' |tr '\t' '\n' |less -S

4.从gff或者gft文件中获取基因的ID与symbol对应关系，以及biotype类型

方法一：

1. zless -S Homo_sapiens.GRCh38.113.chr.gtf.gz: 使用 zless 查看压缩的 GTF 文件 Homo_sapiens.GRCh38.113.chr.gtf.gz 内容。-S 表示启用水平滚动，以便长行可以在屏幕上滚动查看。
2. awk -F'\t' '{if(3=="gene"){print$9}}': 使用 awk 来处理每一行，-F'\t' 表示以制表符（\t）作为分隔符。3=="gene"：仅选择第 3 列等于 "gene" 的行，这意味着我们只提取类型为基因的数据。print $9：打印第 9 列，这一列通常包含基因的详细注释信息，例如 gene_id 和 gene_name。
3. awk -F';' '{print 3,5}': 进一步用 awk 处理数据，-F';' 表示以分号（;）作为分隔符。print 1,5：选择第 1、3 和 5 列，这些列通常包含 gene_id 和 gene_name 等信息。
4. awk '{print 4"\t"$6}': 继续用 awk 对之前的输出进行处理。打印第 2、4 和 6 列，并在它们之间用 \t 制表符分隔，提取所需的字段。
5. sed 's/"//g': 使用 sed 删除输出中的所有双引号（"），s/"//g 表示将双引号替换为空字符。
6. grep 'protein_coding': 使用 grep 筛选仅包含 "protein_coding" 的行，这样只保留编码蛋白质的基因信息。
7. protein_coding_id2name.txt: 将最终结果重定向并保存到 protein_coding_id2name.txt 文件中。

方法二跟方法一是类似的。

方法三：

1. zless Homo_sapiens.GRCh38.113.chr.gtf.gz: 使用 zless 查看压缩的 GTF 文件内容，行为与前面相同。
2. grep -v "#": 使用 grep 过滤掉以 # 开头的注释行，-v 表示反选，即只保留非注释行。
3. awk -F '\t' '{if($3=="gene") {print $9}}': 使用 awk 处理以制表符分隔的字段，只选择第 3 列为 "gene" 的行，并打印第 9 列（基因信息）。
4. cut -d";" -f1,3,5: 使用 cut 命令，以 ; 作为分隔符，选择第 1、3 和 5 字段，类似于之前的 awk -F';' 操作。
5. cut -d' ' -f2,4,6: 再次使用 cut 命令，以空格作为分隔符，选择第 2、4 和 6 字段，这样就提取了所需信息。
6. sed 's#"##g': 使用 sed 去掉双引号，和之前的 sed 命令作用相同。
7. grep "protein_coding": 筛选编码蛋白质的基因信息。
8. protein_coding_id2name.txt: 将结果保存到 protein_coding_id2name.txt 文件中。

# 从gff或者gft文件中获取ID与symbol对应关系，以及biotype类型
# gtf=/Desktop/RNA/Human-3-NPC-Tra/Homo_sapiens.GRCh38.113.gtf.gz
zless -S Homo_sapiens.GRCh38.113.gtf.gz | awk -F'\t' '{if($3=="gene"){print$9}}' |awk -F';' '{print$1,$3,$5}' |awk '{print$2"\t"$4"\t"$6}' |sed 's/"//g' |grep 'protein_coding' > protein_coding_id2name.txt

zless -S $gtf| awk -F'\t' '{if($3=="gene"){print$9}}' |awk -F';' '{print$1,$3,$5}' |awk '{print$2"\t"$4"\t"$6}' |sed 's/"//g' |grep 'protein_coding' > protein_coding_id2name.txt

# 或者
zless Homo_sapiens.GRCh38.113.chr.gtf.gz | grep -v "#" | awk -F '\t' '{if($3=="gene") {print $9}}' | cut -d";" -f1,3,5| cut -d' ' -f2,4,6 | sed 's#"##g' | grep "protein_coding" > protein_coding_id2name.txt

5.Hisat2比对

HISAT2（Hierarchical Indexing for Spliced Alignment of Transcripts 2）是一个用于比对 RNA-Seq 数据到参考基因组的高效和快速的比对工具。它是 HISAT 的升级版本，专门设计来处理高通量测序产生的大规模数据集，特别是 RNA-Seq 数据。

使用 hisat2 进行比对，需要对参考基因组进行构建索引。

# cd ~/database/genome
# 其中 fa 文件是参考基因组，后面的是生成索引文件的前缀
# ch是人的小鼠是m
# 下载的时候是gz文件，需要把gz文件解压缩
hisat2-build  Homo_sapiens.GRCh38.dna.primary_assembly.fa  GRCh38.dna

# 构建好索引，会生成多个 ht 文件
$ ls
GRCh38.dna.1.ht2
GRCh38.dna.2.ht2
GRCh38.dna.3.ht2
GRCh38.dna.4.ht2
GRCh38.dna.5.ht2
GRCh38.dna.6.ht2
GRCh38.dna.7.ht2
GRCh38.dna.8.ht2

Homo_sapiens.GRCh38.dna.primary_assembly.fa

正式比对命令：

里面的路径信息需要更换为使用者自己的。

创建好文件夹，把中间文件放进去

1. -p 2：指定使用两个线程进行并行计算，以提高处理速度。
2. -x ：选项指定了参考基因组的索引文件的前缀。{index} 已经生成的索引文件的路径（例如，GRCh38.dna）。
3. -1 ./trim_galore/SRR23881762_1_val_1.fq 和 -2 ./trim_galore/SRR23881762_2_val_2.fq：这些选项指定了两个输入文件，这里是配对末端序列数据（paired-end reads），分别是 1 和 2 端的数据。文件经过 trim_galore 处理过，因此名称包含 _val_1 和 _val_2 后缀。
4. |：管道符，表示将 hisat2 的输出直接传递给下一个命令，即 samtools，不生成中间文件。
5. samtools sort -@ 2 -o ./hisat2/SRR23881762.Hisat_aln.sorted.bam -：samtools sort：用于对 BAM 文件（对齐结果的二进制格式）进行排序。-@ 2：指定使用两个线程进行排序。-o ./hisat2/SRR23881762.Hisat_aln.sorted.bam：指定输出排序后的 BAM 文件，并保存到 ./hisat2 目录中，文件名为 SRR23881762.Hisat_aln.sorted.bam。-：表示从管道中读取数据。

# 返回项目工作目录
cd ~/Desktop/RNA/Human-3-NPC-Tra/

# 赋值索引前缀(就是生成的8个文件)
index=~/Desktop/RNA/Human-3-NPC-Tra/GRCh38.113/hisat2Index/GRCh38.dna

# 单个样本比对
#hisat2 -p 2  -x  ${index}                 \
#-1 ./trim_galore/SRR23881762_1_val_1.fq \
#-2 ./trim_galore/SRR23881762_2_val_2.fq \
#-S ./hisat2/SRR23881762.Hisat_aln.sam

# sam 转 bam，并且进行 sort
#samtools sort -@ 2  ./hisat2/SRR23881762.Hisat_aln.sam -o ./hisat2/SRR23881762.Hisat_aln.sorted.bam

# 通常来说，不需要生成 sam 文件的，上面几句代码可以通过管道符 | 合并为一句，最后需要有占位符 - 
hisat2 -p 2  -x  ${index}                 \
-1 ./trim_galore/SRR23881762_1_val_1.fq \
-2 ./trim_galore/SRR23881762_2_val_2.fq \
 | samtools sort -@ 2 -o ./hisat2/SRR23881762.Hisat_aln.sorted.bam - 

# 对bam建索引
# 这个命令会生成一个 .bai 索引文件，这个索引文件保存了数据在 BAM 文件中的位置信息，以便快速查询。
samtools index ./hisat2/SRR23881762.Hisat_aln.sorted.bam ./hisat2/SRR23881762.Hisat_aln.sorted.bam.bai

得到一个比对结果

6.多样本，把下面代码保存为 Hisat.sh 脚本

index=~/Desktop/RNA/Human-3-NPC-Tra/GRCh38.113/Hisat2Index/GRCh38.dna

cat sample.ID | while read id
do
  hisat2 -p 2 -x ${index}              \
  -1 ./trim_galore/${id}_1_val_1.fq.gz \
  -2 ./trim_galore/${id}_2_val_2.fq.gz \
  2>${id}.log  |                       \
  samtools sort -@ 2 -o ./hisat2/${id}.Hisat_aln.sorted.bam -  
  
done

7.提交后台运行

nohup  bash  Hisat.sh  >Hisat.log  &

数据表达定量

featureCounts 是一个用于对对齐的 RNA-Seq 数据进行特征计数的工具，常用于生成基因或转录本表达矩阵。

1. -a ：选项指定注释文件，这里使用了之前定义的文件路径gtf。
2. -o ./featureCounts/all.count.txt：-o 选项指定输出文件，这里输出的是 ./featureCounts/all.count.txt。
3. -p：用于配对末端数据（paired-end reads）的计数。
4. -T 6：指定使用 6 个线程以提高处理速度。
5. -t exon：-t 选项指定要计数的特征类型，这里是 exon（外显子）。
6. -g gene_id：-g 选项指定要使用的分组特征，这里是 gene_id（基因ID）。
7. ./hisat2/*.sorted.bam：指定输入的 BAM 文件，这里使用的是 ./hisat2 目录下所有 .sorted.bam 结尾的文件。
8. grep -v '#' ./featureCounts/all.count.txt：过滤掉 all.count.txt 文件中的注释行，这些行通常以 # 开头。-v 选项表示反转匹配（即排除以 # 开头的行）。
9. cut -f 1,7-：cut 命令用于提取特定列，这里提取的是第 1 列（通常是基因 ID）和第 7 列开始的所有列（通常是样本的计数数据）。
10. sed "s@./hisat2/@@g"：使用 sed 替换文本，s@./hisat2/@@g 表示将路径 ./hisat2/ 替换为空字符串（即删除它）。@ 是分隔符，可以用其他符号代替。
11. sed 's#.Hisat_aln.sorted.bam##g'：将 BAM 文件的后缀 .Hisat_aln.sorted.bam 替换为空字符串，这样可以得到干净的样本名称。
12. ./featureCounts/raw_counts.txt：将结果重定向到 ./featureCounts/raw_counts.txt 文件。
13. head ./featureCounts/raw_counts.txt：查看 raw_counts.txt 文件的前 10 行。
14. | column -t：column 命令将输出的数据进行格式化对齐显示，-t 选项会使用空格对列进行对齐，以提高可读性。

# 返回项目工作目录
cd ~/Desktop/RNA/Human-3-NPC-Tra/

# featureCounts对bam文件进行计数

## 定义输入输出文件
gtf=~/Desktop/RNA/Human-3-NPC-Tra/GRCh38.113/Homo_sapiens.GRCh38.113.chr.gtf.gz

# 如果是Version v2.0.3，使用下面的代码：
featureCounts -a $gtf -o ./featureCounts/all.count.txt -p -T 6 -t exon -g gene_id ./hisat2/*.sorted.bam

# 得到表达矩阵
# 处理表头，要换成自己的路径
grep -v '#' ./featureCounts/all.count.txt |cut -f 1,7- | sed "s@./hisat2/@@g" | sed 's#.Hisat_aln.sorted.bam##g' > ./featureCounts/raw_counts.txt

sed s///
sed s###
sed s%%%


# 列对齐显示
head ./featureCounts/raw_counts.txt  | column -t

既往推文

转录组上游分析流程(一）：https://mp.weixin.qq.com/s/bwUFJ-kBdUTp9WyQRQPnyg

转录组上游分析流程(二）：https://mp.weixin.qq.com/s/Lq9lES4M6ZzzleXfg0iWCw

转录组上游分析流程(三）：https://mp.weixin.qq.com/s/EIoPqIKC4IuOFChCmHlsVQ

参考资料：

1. Ensembl官网：https://asia.ensembl.org/index.html
2. 生信技能树：https://mp.weixin.qq.com/s/dxoorMYHU-tlxMcICnTQTg
3. 生信菜鸟团：https://mp.weixin.qq.com/s/Z5YNRkJ6tlY_tAeuUfL8SA

致谢：感谢曾老师/新叶老师以及生信技能树团队全体成员。

注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟

- END -