c
chang
V1
2023/01/09阅读：316主题：默认主题
RNA-seq上游数据处理
软件安装
此次需要用到的有
-
质控 fastqc
-
质控过滤 trim-galore
使用
trim_galore命令 -
RNA序列对比软件 hisat2
-
序列比对 samtools
-
计算表达量（定量工具）subread
这里使用的是
featureCounts -
转录本定量 salmon
安装代码如下
conda create -n rna
conda activate rna
conda install -y -c bioconda fastqc trim-galore hisat2 subread
conda install -y -c bioconda salmon # salmon-0.14.2
conda install -y -c bioconda samtools # samtools-1.6
# 如果是 STAR-Fusion 看融合基因，RNA编辑位点
# 可能性需要为 STAR-Fusion 专门的设置一个conda小环境，它毕竟特殊。
验证代码如下
fastqc --help 1>/dev/null
trim_galore --help 1>/dev/null ###注意：这里是下划线 不是横线
hisat2 --help 1>/dev/null
featureCounts --help 1>/dev/null
下载样本数据
使用prefetch下载
#下载单个数据
prefetch SRR11618610
#如果同时下载两个或多个，并行即可
prefetch SRR11618610 SRR11618616
#批量全部下载，样本量大时，需谨慎服务器资源
cat SRR_Acc_List.txt |while read id;do (prefetch -X 100G $id &);done
sra数据转化成fastq格式
注意路径的使用
-O输出路径 -e线程数 -A输入文件
*表示任何长度和数值的通配符 “TEXTGROUP”文件夹为测试的输出文件夹
-
批处理方法
ls -d SRR/*sra | while read id;do ( fasterq-dump -O ./TEXTGROUP --split-files -e 2 ./$id --include-technical & );done -
单独处理的方法
*.sra-> *.fastq
fastq-dump --split-3 ~/SSR/SRR11618621.sra
fastq-dump --split-3 ~/SSR/SRR11618616.sra
(download) tree116@st8-Horsea-12U:~$ fastq-dump --split-3 -O ~/text -A ~/SSR/*.sra
Read 676511 spots for /home/st8/ssd2/tree116/SSR/SRR11618616.sra
Written 676511 spots for /home/st8/ssd2/tree116/SSR/SRR11618616.sra
Read 832516 spots for /home/st8/ssd2/tree116/SSR/SRR11618621.sra
Written 832516 spots for /home/st8/ssd2/tree116/SSR/SRR11618621.sra
(download) tree116@st8-Horsea-12U:~$ ls text
_home_st8_ssd2_tree116_SSR_SRR11618610.sra_1.fastq
_home_st8_ssd2_tree116_SSR_SRR11618610.sra_2.fastq
将fastq文件压缩
ls TEXTGROUP/*fastq |while read id;do ( gzip $id &);done
(rna) [tree116@st8-Horsea-12U 17:01:31 ~ ]$ ls TEXTGROUP/*fastq |while read id;do ( gzip $id &);done
(rna) [tree116@st8-Horsea-12U 17:01:48 ~ ]$ ls TEXTGROUP
SRR11618610_1.fastq.gz SRR11618616_1.fastq.gz SRR11618621_1.fastq.gz
SRR11618610_2.fastq SRR11618616_2.fastq SRR11618621_2.fastq
SRR11618610_2.fastq.gz SRR11618616_2.fastq.gz SRR11618621_2.fastq.gz
SRR11618616_1.fastq SRR11618621_1.fastq
质控
-
使用fastqc进行质控
nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log & -
使用multiqc对fastqc结果进行整合
nohup multiqc ./*.zip -o ./multiqc > ./multiqc.log & -
trim_galore对测序结果进行处理（但是运行这个命令时遇到了麻烦）
ls *gz |while read id;do (nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o ./ $id & );done
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Multicore support not enabled. Proceeding with single-core trimming.
Path to Cutadapt set as: 'cutadapt' (default)
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt seems to be working fine (tested command 'cutadapt --version')
Cutadapt version: 4.1
single-core operation.
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
Cutadapt version: 4.1
single-core operation.
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
Cutadapt version: 4.1
single-core operation.
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
Cutadapt version: 4.1
single-core operation.
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
Cutadapt version: 4.1
single-core operation.
Cutadapt version: 4.1
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
single-core operation.
Please provide an even number of input files for paired-end FastQ trimming! Aborting ...
#觉得是双端变单端，输入参数不正确
将--paired参数去掉，测试单端比对，这次出了结果，说明单次输入文件应为两个
(rna) [tree116@st8-Horsea-12U 11:30:42 ~/SRR/test ]$ ls out -l
total 141468
-rw-rw-r-- 1 tree116 tree116 3002 12月 5 11:30 SRR11618610_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6371450 12月 5 11:30 SRR11618610_1_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3026 12月 5 11:30 SRR11618610_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6758028 12月 5 11:30 SRR11618610_2_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3443 12月 5 11:30 SRR11618616_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 29130386 12月 5 11:30 SRR11618616_1_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3432 12月 5 11:30 SRR11618616_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 30419370 12月 5 11:30 SRR11618616_2_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3043 12月 5 11:30 SRR11618621_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 35856265 12月 5 11:30 SRR11618621_1_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3094 12月 5 11:30 SRR11618621_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 36291370 12月 5 11:30 SRR11618621_2_trimmed.fq.gz
(rna) [tree116@st8-Horsea-12U 11:55:42 ~/SRR/test/out ]$ cat SRR11618610_1.fastq.gz_trimming_report.txt
SUMMARISING RUN PARAMETERS
==========================
Input filename: SRR11618610_1.fastq.gz
Trimming mode: paired-end
Trim Galore version: 0.6.7
Cutadapt version: 4.1
Number of cores used for trimming: 1
Quality Phred score cutoff: 25
Quality encoding type selected: ASCII+33
Using Nextera adapter for trimming (count: 452). Second best hit was smallRNA (count: 2)
Adapter sequence: 'CTGTCTCTTATA' (Nextera Transposase sequence; auto-detected)
Maximum trimming error rate: 0.1 (default)
Minimum required adapter overlap (stringency): 3 bp
Minimum required sequence length for both reads before a sequence pair gets removed: 36 bp
Output file will be GZIP compressed
但是样本是双端，修改了批处理命令，参考：https://cloud.tencent.com/developer/article/1703054
for i in SRR11618610 SRR11618616 SRR11618621; do trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired /home/st8/ssd2/tree116/fastqc/${i}_1.fastq.gz /home/st8/ssd2/tree116/fastqc/${i}_2.fastq.gz --gzip -o /home/st8/ssd2/tree116/fastqc/trim_galoreR/; done
（ 12月6日更新，群内tree011小伙伴分享了一个更简单的命令
ls *gz |xargs -n 2 -t nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -j 4 -o ./ &
管道符后面的xargs指定传两个文件名作为参数，默认"_1"和"_2"是排好序的，可以直接使用。-t参数可以显示即将执行的代码 ）
在文件夹内成功输出
(rna) [tree116@st8-Horsea-12U 18:39:09 ~/fastqc ]$ for i in SRR11618610 SRR11618616 SRR11618621; do trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired /home/st8/ssd2/tree116/fastqc/${i}_1.fastq.gz /home/st8/ssd2/tree116/fastqc/${i}_2.fastq.gz --gzip -o /home/st8/ssd2/tree116/fastqc/trim_galoreR/; done
(rna) [tree116@st8-Horsea-12U 18:39:00 ~/fastqc ]$ ls trim_galoreR -ll
total 140612
-rw-rw-r-- 1 tree116 tree116 3022 12月 5 18:37 SRR11618610_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6293855 12月 5 18:37 SRR11618610_1_val_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 3245 12月 5 18:37 SRR11618610_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6717968 12月 5 18:37 SRR11618610_2_val_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 3462 12月 5 18:37 SRR11618616_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 28780467 12月 5 18:38 SRR11618616_1_val_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 3652 12月 5 18:38 SRR11618616_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 30172038 12月 5 18:38 SRR11618616_2_val_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 3062 12月 5 18:38 SRR11618621_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 35755068 12月 5 18:39 SRR11618621_1_val_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 3313 12月 5 18:39 SRR11618621_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 36224861 12月 5 18:39 SRR11618621_2_val_2.fq.gz
这里看到了一种用脚本处理的方法 尝试一下 https://www.jianshu.com/p/74a328f4d23a
grep _1.fastq.gz>gz1
grep _2.fastq.gz>gz2
#创建名为config的矩阵，以1，2为两列，以进行分组
paste gz1 gz2>config
(rna) [tree116@st8-Horsea-12U 14:52:10 ~/SRR/test ]$ cat config
SRR11618610_1.fastq.gz SRR11618610_2.fastq.gz
SRR11618616_1.fastq.gz SRR11618616_2.fastq.gz
SRR11618621_1.fastq.gz SRR11618621_2.fastq.gz
ls /bin/sh -al结果bin/sh指向为dash 而我们的脚本为bash脚本，所以运行时应使用bash命令运行
#创建脚本文件
vi trim.sh
#以下为.sh文件内容
#!/bin/bash
#运行时也需要使用bash，否则会出现`Bad substitution`的问题
dir=/home/st8/ssd2/tree116/SRR/test/outsh/
cat config |while read id
do
arr=${id}
fq1=${arr[0]}
fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired $fq1 $fq2 -o $dir & > $dir/bashtest.log
done
#运行
bash trim.sh
结果如下
(rna) [tree116@st8-Horsea-12U 16:22:56 ~/SRR/test/outsh ]$ ls -alh
total 94M
drwxrwxrwx 2 tree116 tree116 4.0K 12月 5 16:23 .
drwxrwxr-- 6 tree116 tree116 4.0K 12月 5 16:15 ..
-rw-rw-r-- 1 tree116 tree116 0 12月 5 16:22 bashtest.log
-rw-rw-r-- 1 tree116 tree116 2.9K 12月 5 16:22 SRR11618610_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6.1M 12月 5 16:22 SRR11618610_1_val_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 3.1K 12月 5 16:22 SRR11618610_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6.5M 12月 5 16:22 SRR11618610_2_val_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 3.3K 12月 5 16:22 SRR11618616_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 28M 12月 5 16:22 SRR11618616_1_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 670 12月 5 16:22 SRR11618616_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 13M 12月 5 16:23 SRR11618616_2_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 3.0K 12月 5 16:23 SRR11618621_1.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 35M 12月 5 16:23 SRR11618621_1_trimmed.fq.gz
-rw-rw-r-- 1 tree116 tree116 670 12月 5 16:23 SRR11618621_2.fastq.gz_trimming_report.txt
-rw-rw-r-- 1 tree116 tree116 6.0M 12月 5 16:23 SRR11618621_2_trimmed.fq.gz
-
使用fastqc重新质控
nohup fastqc -t 6 -o ./ SRR*.fq.gz >qc.log &这里的格式为fq.gz，上文为fastq.gz
nohup multiqc ./*.zip -o ./multiqc > ./multiqc.log &
质控'
fastp直接质控(双端) 参考： https://zhuanlan.zhihu.com/p/33601691
for i in SRR11618610 SRR11618616 SRR11618621
do
nohup fastp -i ../fastqc/${i}_1.fastq.gz -o ${i}_1.fq.gz -I ../fastqc/${i}_2.fastq.gz -O ${i}_2.fq.gz &
done
(rna) [tree116@st8-Horsea-12U 20:42:42 ~/RNA/fastp ]$ ls -lh
total 1.2G
-rw-rw-r-- 1 tree116 tree116 448K 12月 30 20:42 fastp.html
-rw-rw-r-- 1 tree116 tree116 111K 12月 30 20:42 fastp.json
-rw------- 1 tree116 tree116 3.1K 12月 30 20:42 nohup.out
-rw-rw-r-- 1 tree116 tree116 35M 12月 5 21:10 SRR11618610_1.fq
-rw-rw-r-- 1 tree116 tree116 6.9M 12月 30 20:42 SRR11618610_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 35M 12月 5 21:10 SRR11618610_2.fq
-rw-rw-r-- 1 tree116 tree116 7.2M 12月 30 20:42 SRR11618610_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 169M 12月 5 21:10 SRR11618616_1.fq
-rw-rw-r-- 1 tree116 tree116 32M 12月 30 20:42 SRR11618616_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 169M 12月 5 21:10 SRR11618616_2.fq
-rw-rw-r-- 1 tree116 tree116 33M 12月 30 20:42 SRR11618616_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 218M 12月 5 21:10 SRR11618621_1.fq
-rw-rw-r-- 1 tree116 tree116 41M 12月 30 20:42 SRR11618621_1.fq.gz
-rw-rw-r-- 1 tree116 tree116 218M 12月 5 21:10 SRR11618621_2.fq
-rw-rw-r-- 1 tree116 tree116 41M 12月 30 20:42 SRR11618621_2.fq.gz
-rw-rw-r-- 1 tree116 tree116 218M 12月 5 22:11 SRR11618621_2.fq.gzfor
输出文件nohup.out如下：
Read1 before filtering:
total reads: 137124
total bases: 13712400
Q20 bases: 13189766(96.1886%)
Q30 bases: 12571273(91.6781%)
Read2 before filtering:
total reads: 137124
total bases: 13712400
Q20 bases: 12892522(94.0209%)
Q30 bases: 12030598(87.7352%)
Read1 after filtering:
total reads: 131840
total bases: 13147049
Q20 bases: 12746851(96.956%)
Q30 bases: 12183232(92.6689%)
Read2 after filtering:
total reads: 131840
total bases: 13147049
Q20 bases: 12522143(95.2468%)
Q30 bases: 11724352(89.1786%)
Filtering result:
reads passed filter: 263680
reads failed due to low quality: 10552
reads failed due to too many N: 16
reads failed due to too short: 0
reads with adapter trimmed: 2116
bases trimmed due to adapters: 74124
Duplication rate: 35.6575%
Insert size peak (evaluated by paired-end reads): 129
JSON report: fastp.json
HTML report: fastp.html
fastp -i ../fastqc/SRR11618610_1.fastq.gz -o SRR11618610_1.fq.gz -I ../fastqc/SRR11618610_2.fastq.gz -O SRR11618610_2.fq.gz
fastp v0.23.2, time used: 1 seconds
Read1 before filtering:
total reads: 676511
total bases: 67651100
Q20 bases: 64491682(95.3298%)
Q30 bases: 61523379(90.9422%)
Read2 before filtering:
total reads: 676511
total bases: 67651100
Q20 bases: 63554856(93.945%)
Q30 bases: 59823489(88.4294%)
Read1 after filtering:
total reads: 640856
total bases: 63396915
Q20 bases: 61359387(96.7861%)
Q30 bases: 58807105(92.7602%)
Read2 after filtering:
total reads: 640856
total bases: 63396915
Q20 bases: 60742378(95.8128%)
Q30 bases: 57477207(90.6625%)
Filtering result:
reads passed filter: 1281712
reads failed due to low quality: 71226
reads failed due to too many N: 84
reads failed due to too short: 0
reads with adapter trimmed: 42276
bases trimmed due to adapters: 1384228
Duplication rate: 55.2044%
Insert size peak (evaluated by paired-end reads): 129
JSON report: fastp.json
HTML report: fastp.html
fastp -i ../fastqc/SRR11618616_1.fastq.gz -o SRR11618616_1.fq.gz -I ../fastqc/SRR11618616_2.fastq.gz -O SRR11618616_2.fq.gz
fastp v0.23.2, time used: 3 seconds
Read1 before filtering:
total reads: 832516
total bases: 83251600
Q20 bases: 81395406(97.7704%)
Q30 bases: 78398402(94.1704%)
Read2 before filtering:
total reads: 832516
total bases: 83251600
Q20 bases: 81146923(97.4719%)
Q30 bases: 78000409(93.6924%)
Read1 after filtering:
total reads: 824867
total bases: 82295395
Q20 bases: 80606696(97.948%)
Q30 bases: 77688078(94.4015%)
Read2 after filtering:
total reads: 824867
total bases: 82295395
Q20 bases: 80454675(97.7633%)
Q30 bases: 77395871(94.0464%)
Filtering result:
reads passed filter: 1649734
reads failed due to low quality: 15190
reads failed due to too many N: 108
reads failed due to too short: 0
reads with adapter trimmed: 14158
bases trimmed due to adapters: 382916
Duplication rate: 42.696%
Insert size peak (evaluated by paired-end reads): 169
JSON report: fastp.json
HTML report: fastp.html
fastp -i ../fastqc/SRR11618621_1.fastq.gz -o SRR11618621_1.fq.gz -I ../fastqc/SRR11618621_2.fastq.gz -O SRR11618621_2.fq.gz
fastp v0.23.2, time used: 3 seconds
比对
首先是人类参考基因组和注释文件 在hisat2官网下载 https://daehwankimlab.github.io/hisat2/download/#h-sapiens
这里使用Ch38中的 genome_tran 
下载后用FTP拖进服务器
使用hisat2进行比对
for i in 10 16 21
do
hisat2 --dta -t -x ./grch38/genome -1 ./RNA/fastp/SRR116186${i}_1.fq.gz -2 RNA/fastp/SRR116186${i}_2.fq.gz -S ./RNA/hisat2/SRR116186${i}.sam &
done
#结果如下
(rna) [tree116@st8-Horsea-12U 17:13:52 ~/RNA/h2test ]$ ls -lh
total 1.2G
-rw-rw-r-- 1 tree116 tree116 90M 12月 29 17:11 SRR11618610.sam
-rw-rw-r-- 1 tree116 tree116 481M 12月 29 17:11 SRR11618616.sam
-rw-rw-r-- 1 tree116 tree116 616M 12月 29 17:11 SRR11618621.sam
输出如下 overall alignment rate分别为86.48% 66.24% 82.24%
Time loading forward index: 00:00:02
Time loading forward index: 00:00:02
Time loading forward index: 00:00:02
Time loading reference: 00:00:01
Time loading reference: 00:00:01
Time loading reference: 00:00:01
Multiseed full-index search: 00:00:09
131840 reads; of these:
131840 (100.00%) were paired; of these:
31540 (23.92%) aligned concordantly 0 times
97199 (73.72%) aligned concordantly exactly 1 time
3101 (2.35%) aligned concordantly >1 times
----
31540 pairs aligned concordantly 0 times; of these:
2189 (6.94%) aligned discordantly 1 time
----
29351 pairs aligned 0 times concordantly or discordantly; of these:
58702 mates make up the pairs; of these:
35646 (60.72%) aligned 0 times
20843 (35.51%) aligned exactly 1 time
2213 (3.77%) aligned >1 times
86.48% overall alignment rate
Time searching: 00:00:11
Overall time: 00:00:13
Multiseed full-index search: 00:00:44
640856 reads; of these:
640856 (100.00%) were paired; of these:
254355 (39.69%) aligned concordantly 0 times
344661 (53.78%) aligned concordantly exactly 1 time
41840 (6.53%) aligned concordantly >1 times
----
254355 pairs aligned concordantly 0 times; of these:
3482 (1.37%) aligned discordantly 1 time
----
250873 pairs aligned 0 times concordantly or discordantly; of these:
501746 mates make up the pairs; of these:
432723 (86.24%) aligned 0 times
52676 (10.50%) aligned exactly 1 time
16347 (3.26%) aligned >1 times
66.24% overall alignment rate
Time searching: 00:00:45
Overall time: 00:00:47
Multiseed full-index search: 00:00:53
824867 reads; of these:
824867 (100.00%) were paired; of these:
191323 (23.19%) aligned concordantly 0 times
587332 (71.20%) aligned concordantly exactly 1 time
46212 (5.60%) aligned concordantly >1 times
----
191323 pairs aligned concordantly 0 times; of these:
4757 (2.49%) aligned discordantly 1 time
----
186566 pairs aligned 0 times concordantly or discordantly; of these:
373132 mates make up the pairs; of these:
293021 (78.53%) aligned 0 times
66181 (17.74%) aligned exactly 1 time
13930 (3.73%) aligned >1 times
82.24% overall alignment rate
Time searching: 00:00:54
Overall time: 00:00:56
sam 转化成bam
#转换
for i in 10 16 21
do
samtools view SRR116186${i}.sam -b > SRR116186${i}.bam #转换
samtools sort SRR116186${i}.bam -o SRR116186${i}_sorted.bam #排序
samtools index SRR116186${i}_sorted.bam; #索引
done
查看单个结果比对后的总体质量 给出BAM文件的比对结果
samtools flagstat SRR11618610_sorted.bam
(rna) [tree116@st8-Horsea-12U 21:35:48 ~/RNA/fastp ]$ samtools flagstat SRR11618610_sorted.bam
288817 + 0 in total (QC-passed reads + QC-failed reads)
25137 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
253171 + 0 mapped (87.66% : N/A)
263680 + 0 paired in sequencing
131840 + 0 read1
131840 + 0 read2
200600 + 0 properly paired (76.08% : N/A)
207422 + 0 with itself and mate mapped
20612 + 0 singletons (7.82% : N/A)
2340 + 0 with mate mapped to a different chr
1979 + 0 with mate mapped to a different chr (mapQ>=5)
(rna) [tree116@st8-Horsea-12U 21:38:51 ~/RNA/fastp ]$ samtools flagstat SRR11618616_sorted.bam
1610889 + 0 in total (QC-passed reads + QC-failed reads)
329177 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
1178166 + 0 mapped (73.14% : N/A)
1281712 + 0 paired in sequencing
640856 + 0 read1
640856 + 0 read2
773002 + 0 properly paired (60.31% : N/A)
784814 + 0 with itself and mate mapped
64175 + 0 singletons (5.01% : N/A)
6486 + 0 with mate mapped to a different chr
5589 + 0 with mate mapped to a different chr (mapQ>=5)
(rna) [tree116@st8-Horsea-12U 21:43:24 ~/RNA/fastp ]$ samtools flagstat SRR11618621_sorted.bam
1983306 + 0 in total (QC-passed reads + QC-failed reads)
333572 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
1690285 + 0 mapped (85.23% : N/A)
1649734 + 0 paired in sequencing
824867 + 0 read1
824867 + 0 read2
1267088 + 0 properly paired (76.81% : N/A)
1285356 + 0 with itself and mate mapped
71357 + 0 singletons (4.33% : N/A)
6374 + 0 with mate mapped to a different chr
4009 + 0 with mate mapped to a different chr (mapQ>=5)
(rna) [tree116@st8-Horsea-12U 21:43:32 ~/RNA/fastp ]$
samtools sort -n test.bam -o sort_left.bam -n表示左
samtools view sort_left.bam | head
samtools view SRR11618621_sorted.bam 1:1234-123456 | head
对BAM文件进行统计分析
①RSeQC 三个文件的统计分析如下图：
bam_stat.py -i SRR11618610_sorted.bam
②使用bed文件作为参考 查看基因组覆盖率
需要提供bed文件，可以直接RSeQC的网站下载（本文统一使用hg38），或者可以用gtf转换 https://sourceforge.net/projects/rseqc/postdownload
2023.1.1日凌晨更新：事实上，此时定量已经成功，且对比率水平不错(如下所示，具体定量细节见下文“ 定量 ”部分)，怀疑是bed文件出了问题。（未验证）
Load annotation file gencode.v42.annotation.gtfFeatures : 1643237
Meta-features : 62696
Chromosomes/contigs : 25
Process BAM file SRR11618610 sorted.bam...
Paired-end reads are included
Total alignments : 147824
Successfully assigned alignments : 57279 (38.7%)
Running time : 0.01 minutes
Process BAM file SRR11618616 sorted.bam...
Paired-end reads are included.
Total alignments : 831973
Successfully assigned alignments : 308850 (37.1%)
Running time : 0.01 minutes
Process BAM file SRR11618621 sorted.bam...
Paired-end reads are included.
Total alignments : 1012896
Successfully assigned alignments : 485746 (48.0%)Running time : 0.01 minutes
定量
这里使用featureCounts工具 https://www.bioinfo-scrounger.com/archives/407/
gtf相关阅读： https://mp.weixin.qq.com/s/YHWLcZYeKLEMufUS-TLHVQ
代码举例如下，这里主要讨论不同gtf带来的表达矩阵的结果差异及其原因。
gtf=$HOME/database/GRCh38.106/Homo_sapiens.GRCh38.106.chr.gtf
nohup featureCounts -T 5 -p -t exon -g gene_id -a $gtf \
-o all.id.txt *sorted.bam 1>counts.id.log 2>&1 &
下载参考gtf文件的方法
https://blog.csdn.net/u011262253/article/details/89363809 https://mp.weixin.qq.com/s/O9KZdU9XvqW0_ZWtJM-rnA
各个参数的说明：
-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对paired-end数据
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数
-
首先NCBI的gff文件 GRCh38_latest_genomic.gff
nohup featureCounts -T 5 -p -t exon -g ID -a $gtf -o resulNCBI/all.id.txt *sorted.bam 1> resulNCBI/counts.id.log 2>&1 &
-g 当提供参考的gtf/gff的时候，我们需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计，默认为“gene_id”，注意！选择gtf/gff中提供的id identifier！！！GFF文件来自NCBI官网，该参数值改为“ID”(见GTF/GFF文件第9列（参考群内同学文章） 结果如下：counts数全部为0，且比对率也为0
-
GENECODE
https://www.gencodegenes.org/human/releases.html gencode.v42.annotation.gtf 
Homo_sapiens.GRCh38.86.chr_patch_hapl_scaff.gtf
使用R查看gtf文件
nohup featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o resultHomo/allHomo.id.txt *sorted.bam 1> resultHomo/countsHomo.id.log 2>&1 & 
后续数据处理参考： https://blog.csdn.net/leo12354/article/details/107536729
这里可以在终端内将需要的几列单独保留下来（包括gene_id 和 各样本表达量） cut -f 1,7,8,9 allHomo.id.txt |grep -v '^#' > feacturCounts.txt
R语言内处理部分：
counts<-read.table("feacturCounts.txt",sep = "\t",header =T)
counts = subset(counts,SRR11618610_sorted.bam != 0 & SRR11618616_sorted.bam != 0 & SRR11618621_sorted.bam != 0) #去掉表达量为0的行
作者介绍
c
chang
V1