samtools view R1.treat.sam | less -S
HWI-ST373:478:C76P5ACXX:1:1101:2241:2250 99 chr2 200203670 27 101M = 200203754 185 TTCTCCATTTACACTTTGTAACTTCACATTCATCCTCTCCATTTACATGGATCATACACACCAAGTAACATCCTCTCCATTTACATAGGGCACATTCCAAG CCCFFFFFHHHHHJJJJJJIJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJIJJJGIJIJJJJJJJEGIJJJJJJJHIJJJJHHHHHHHFFFFDEEEEEEDD AS:i:-11 XS:i:-61 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:38G15C46 YS:i:0 YT:Z:CP
HWI-ST373:478:C76P5ACXX:1:1101:2241:2250 147 chr2 200203754 27 101M = 200203670 -185 ATAGGGCACATTCCAAGTAAATGACTTTGTAACTTCACTTCATTCTCTTTATTTACATAGAGCATACACCAAGTAACCAATGGGAAACCTCTAGAGTATTG DEEEEEEFFFFFFHHHHHHJJJJIJJJJJJJJJJJJIJJJJJJJJJJJIIGIJJIJJJJJJJJJIHHFJJJJJJJJJJJJJJJJJIFJHHHHHFFFFFCCC AS:i:0 XS:i:-52 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:101 YS:i:-11 YT:Z:CP
HWI-ST373:478:C76P5ACXX:1:1101:2282:2106 99 chr7 89102239 42 101M = 89102459 321 GTTATTTCCAAACTGTTTTTTCCATAGAGGTTGAATTAATTTGCATTCCCACCAACAGTATGTGAGCATTCCTTTTTCTTCATATCTGTGCCAATGTTTGT B=BDFFFFHHHHHJIEFIJIJIJIJIJIIJAHIFHIGIIGGGIIJJJGIJJIGIJJIJCHIG@GHIJJJIIJJHHEHHHFFFFFFFFEEEEEEDDDCDDDD AS:i:-4 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:1A99 YS:i:0 YT:Z:CP
HWI-ST373:478:C76P5ACXX:1:1101:2282:2106 147 chr7 89102459 42 101M = 89102239 -321 TGTTTGTTGGTCGCTTGTGTTTCTTCTTTTGAGAAATACCCGTTCATGTCCTTTGCTCAGTTTTGAATGGGGTTGTTTGTTGCATCTTGAGTTGTTTTGGG DDDDDDDDDDDDDDDDDDDDDDEEDEFFFFFFFHHHEAIIJJJJJJJIIIHJIJJJJJJJJJJJJJJJJJJJJJJJJIJJJJJIJJJJHHHHHFFFFFCCC AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:101 YS:i:-4 YT:Z:CP
SAM (Sequence Alignment/Map) 是存储RNA-seq比对结果的格式。BWA、Bowtie/Bowtie2、HISAT/HISAT2、Tophat等均可产生该格式文件。
，以TAB为分割符。SAM格式的设计理念是：灵活存储所有比对信息。SAM一般为11列以上，其中前11列为已定义列，12列及以后为可选列，算附加信息。
由于测序分两种：单端测序(Single-read)和双端测序(Paired-end和Mate-pair)。下文SAM格式中read表示Single-read或Paired-end Read1序列，mate表示Paired-end Read2序列。
image.png
第一列：read name，read的名字通常包括测序平台等信息.
第二列：sum of flags，比对flag数字之和，比对flag用数字表示，不同pdf中对于flag的解释：
image.png
image.png
flag=0：本条read为SE数据，且成功比对到基因组
flag=1：这是PE数据来源的read
flag=2：本条read与本链配对的另外一条read均可以成功比对到参考基因组上
flag=4：本条read不能比对到基因组
flag=8：PE中与本链配对的另外一条read不能比对到基因组
flag=16：本条read是反向互补比对到基因组
flag=32：PE中与本链配对的另外一条read反向互补比对到基因组
flag=64：本条read是R1序列（来自R1.fastq.gz）
flag=128：本条read是R2序列（来自R2.fastq.gz）
flag=256：本条read比对到基因组的多处位置
flag=512：没有通过测序机器本身的质控。这个一般很少见到
flag=1024：PCR重复，猜测判断依据是染色体编号以及起点重点完全一致的reads
flag=2048：存在结构变异，一条read比对到基因组上距离较远的多个位置（可能是不同染色体）
第三列：RNAM，reference sequence name，实际上就是比对到参考序列上的染色体号。若是无法比对，则是*
第四列：position，read比对到参考序列上，第一个碱基所在的位置。若是无法比对，则是0；
第五列：Mapping quality，比对的质量分数，数值越高Reads比对到该位点的可能性越高，255说明此Reads的Mapping quality不可用 (也可简单认为数值越高该read比对到参考基因组上的位置越唯一)；
第六列：CIGAR值，read比对的具体情况，前面的数字代表reads长度
“M”表示 match或 mismatch；
“I”表示 insert；
“D”表示 deletion；
“N”表示 skipped（跳过这段区域）；
“S”表示 soft clipping（被剪切的序列存在于序列中）；
“H”表示 hard clipping（被剪切的序列不存在于序列中）；
“P”表示 padding；
“=”表示 match；
“X”表示 mismatch（错配，位置是一一对应的）；
第七列：MRNM(chr)，mate的reference sequence name，实际上就是mate比对到的染色体号，"="代表与第三列名字一致，若是没有mate，则是*；
第八列：mate position，mate比对到参考序列上的第一个碱基位置，若无mate,则为0；
第九列：ISIZE，Inferred fragment size.详见Illumina中paired end sequencing 和 mate pair sequencing，是负数，推测应该是两条read之间的间隔(待查证)，若无mate则为0；
第十列：Sequence，就是read的碱基序列，如果是比对到互补链上则是reverse completed eg. CGTTTCTGTGGGTGATGGGCCTGAGGGGCGTTCTCN
第十一列：ASCII，read质量的ASCII编码。
第十二列之后：Optional fields，可选的自定义区域：
AS:i 匹配的得分
XS:i 第二好的匹配的得分
YS:i mate 序列匹配的得分
XN:i 在参考序列上模糊碱基的个数
XM:i 错配的个数
XO:i gap open的个数
XG:i gap 延伸的个数
NM:i 经过编辑的序列
YF:i 说明为什么这个序列被过滤的字符串
MD:Z 代表序列和参考序列错配的字符串(例如MD:Z:45A2C3 失配碱基的位点，45,45+2两个位点失配)
XT:A:U read只有一个完整比对，U unique
XT:A:R read有一个以上位置完整比对，R repeat
NM:i:2 read有2个碱基mismatch
X0:i:2 read有2个位置完整比对(与XT:A有对应关系)
X1:i:2 read有2个位置以1个碱基失配比对
XA:Z:Chr3,+1530, 50M,0;Chr4,-7568,50M,1 有0/1个碱基失配的详细比对情况(与XT:A、X0:i、X1:i有对应关系)
image.png
SE数据验证
用简单的数据来验证上面的结论。所以先拿个单端的数据来看。
2.1 统计SE数据常出现的flag
先找一个单端测序的SAM文件来进行分析：
统计这个SAM文件一共都有什么flag存在：
$ samtools view ENCFF000AVM_20220306.sam|awk -F '\t' '{print $2}'|sort|uniq -c
6354014 0
6356967 16
8 20
244687 4
只有4种flag存在。感觉好像有点少，那我们再找一个更大点的文件看看：
统计这个SAM文件一共都有什么flag存在：
$ samtools view ENCFF000AWJ_20220306.sam|awk -F '\t' '{print $2}'|sort|uniq -c
23459767 0
23460459 16
12 20
675438 4
同样还是这4种flag。
看来的单端数据一般只会有这些Flag出现：
0
4
16
20=16+4
flag=0 验证
所用数据包括：
ENCFF000AVM_20220306.sam
ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
我们先找到一条flag=0的read：
$ samtools view ENCFF000AVM_20220306.sam|awk -F '\t' '{if($2==0){print $0}}'|head -1|tr '\t' '\n'
SL-XAU:6:5:1:5:1063:0
0
chr3
15169571
25
36M
*
0
0
GTGCATGCCTGTNGTCTNAGCTACTTAGGACGGTGA
AACCCCCBBCB:%:BB?%7CBBABCBB=BBBB?6B>
XT:A:U
NM:i:2
X0:i:1
X1:i:0
XM:i:2
XO:i:0
XG:i:0
MD:Z:12A4C18
然后去原始的fastq文件中找到原始数据：
$ zgrep -A 2 -B 1 GTGCATGCCTGTNGTCTNAGCTACTTAGGACGGTGA ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
@SL-XAU:6:5:1:5:1063:0/1
GTGCATGCCTGTNGTCTNAGCTACTTAGGACGGTGA
+
AACCCCCBBCB:%:BB?%7CBBABCBB=BBBB?6B>
可以看到原始fastq中的信息和SAM文件中的是一致的。
那么接下来就去参考基因组上拿到SAM文件记录的这段参考基因组序列。
通过CIGAR值我们知道这个read完全比对上去了，那么接下来我们去参考基因组上看看这段序列是什么样子的。首先构建一个bed文件，记录这条read的坐标信息。因为SAM文件坐标是1-based coordinate system，而BED和BAM文件的坐标是0-based coordinate system，所以我们需要将比对的坐标-1。文件 test.bed 内容如下：
chr3 15169570 15169606
利用bedtools工具从参考基因组上进行截取序列：
# 首先需要建立index
$ samtools faidx hg38.fa
# 然后提取目标序列
$ bedtools getfasta -fi hg38.fa -bed test.bed
>chr3:15169570-15169606
gtgcatgcctgtagtctcagctacttaggacggtga
然后这时我们比较下结果：
Reference gtgcatgcctgtagtctcagctacttaggacggtga
sam_read GTGCATGCCTGTNGTCTNAGCTACTTAGGACGGTGA
fq_read GTGCATGCCTGTNGTCTNAGCTACTTAGGACGGTGA
可以看到，这段序列的确是完整比对到了基因组上。
flag=4 验证
我们先找到一条flag=4的read：
$ samtools view ENCFF000AVM_20220306.sam|awk -F '\t' '{if($2==4){print $0}}'|head -1|tr '\t' '\n'
SL-XAU:6:5:1:4:1223:0
4
*
0
0
*
*
0
0
CCAGATGATTCNNNTNTNNTNTTTTTTTTT
BB@ABC?BBC:%%%>%>%%>%:BBBBBCBB
这个很容易就可以看到的确是没有比对上了。不过我们还是找到原始fastq文件中这条read的记录：
$ zgrep -A 2 -B 1 CCAGATGATTCNNNTNTNNTNTTTTTTTTT ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
@SL-XAU:6:5:1:4:1223:0/1
CCAGATGATTCNNNTNTNNTNTTTTTTTTT
+
BB@ABC?BBC:%%%>%>%%>%:BBBBBCBB
flag=16 验证
我们先找到一条flag=16的read：
$ samtools view ENCFF000AVM_20220306.sam|awk -F '\t' '{if($2==16){print $0}}'|head -1|tr '\t' '\n'
SL-XAU:6:5:1:4:1105:0
16
chr19
49699420
25
36M
*
0
0
TGGAGACCAGCCTGGGCANCATANAGACCCCTGTCT
@BBBB<A?ABA>?ABBB;%<AB;%<B<BBBBBBBBA
XT:A:U
NM:i:2
X0:i:1
X1:i:0
XM:i:2
XO:i:0
XG:i:0
MD:Z:18A4G12
我们先用SAM文件中记录的原始序列去找：
$ zgrep -A 2 -B 1 TGGAGACCAGCCTGGGCANCATANAGACCCCTGTCT ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
发现并没有办法找到这条read。接下来我们换成反响互补序列去找：
$ zgrep -A 2 -B 1 AGACAGGGGTCTNTATGNTGCCCAGGCTGGTCTCCA ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
@SL-XAU:6:5:1:4:1105:0/1
AGACAGGGGTCTNTATGNTGCCCAGGCTGGTCTCCA
+
ABBBBBBBB<B<%;BA<%;BBBA?>ABA?A<BBBB@
可以看到，原始fastq文件中序列和SAM文件的结果是反向互补，而原始fastq文件中测序质量和SAM文件的结果是反向的。
那么接下来就去参考基因组上拿到SAM文件记录的这段参考基因组序列。这条read的坐标信息如下：
chr19 49699419 49699455
利用bedtools工具从参考基因组上进行截取序列：
$ bedtools getfasta -fi hg38.fa -bed test.bed
>chr19:49699419-49699455
tggagaccagcctgggcaacatagagacccctgtct
然后这时我们比较下结果：
Reference tggagaccagcctgggcaacatagagacccctgtct
sam_read TGGAGACCAGCCTGGGCANCATANAGACCCCTGTCT
fq_read AGACAGGGGTCTNTATGNTGCCCAGGCTGGTCTCCA
可以看到，原始read（fastq）经过反向互补后变成SAM文件中的read，然后这个read可以完美比对到基因组上chr19:49699419-49699455
flag=20 验证
理论上来说，flag=20=16+4，也就是同时满足：
flag=4：不能比对到基因组
flag=16：反向互补比对到基因组
这就有点矛盾了，在pdf文件中也提到了这种特殊的存在，解释如下：
When 0x4 is set, this indicates whether the unmapped read is stored in its original orientation as it came off the sequencing machine.
还是有点看不懂，所以我们接下来同样的讨论对flag=20的read进行探索：
我们先找到一条flag=20的read：
$ samtools view ENCFF000AVM_20220306.sam|awk -F '\t' '{if($2==20){print $0}}'|head -1|tr '\t' '\n'
SL-XAU:6:5:17:469:1971:0
20
chr19_KI270938v1_alt
1066800
37
36M
*
0
0
GGATCACAGGTCTATCACCCTATTAACCACTCACGG
AB@BA7@BBBBABB@B@==ABBCBCBBCB<@BBCBB
XT:A:U
NM:i:1
X0:i:1
X1:i:0
XM:i:1
XO:i:0
XG:i:0
MD:Z:0T35
因为这个flag包括了16，所以我们用这条read的反向互补序列去获取原始信息：
$ zgrep -A 2 -B 1 CCGTGAGTGGTTAATAGGGTGATAGACCTGTGATCC ENCFF000AVM_20220306_fastp_trimmed.fastq.gz
@SL-XAU:6:5:17:469:1971:0/1
CCGTGAGTGGTTAATAGGGTGATAGACCTGTGATCC
+
BBCBB@<BCBBCBCBBA==@B@BBABBBB@7AB@BA
可以看到，原始fastq文件中序列和SAM文件的结果是反向互补，而原始fastq文件中测序质量和SAM文件的结果是反向的。
那么接下来就去参考基因组上拿到SAM文件记录的这段参考基因组序列。这条read的坐标信息如下：
chr19_KI270938v1_alt 1066799 1066835
利用bedtools工具从参考基因组上进行截取序列：
$ bedtools getfasta -fi hg38.fa -bed test.bed
Feature (chr19_KI270938v1_alt:1066799-1066835) beyond the length of chr19_KI270938v1_alt size (1066800 bp). Skipping.
结果报错了！原因是超过了染色体总的长度1066800
虽然还是对这个flag=20的原因不太清楚，但是可以知道的是，因为它含有flag=4，所以肯定是没法比对上基因组的，可以放心过滤删除！
PE数据验证
我们先找一个双端测序的SAM文件来进行分析：
$ ls -lh 293T_k18ac_rep1.sam
-rw-r--r-- 1 yangliu chenglab 4.1G Mar 22 18:16 293T_k18ac_rep1.bam
统计这个SAM文件一共都有什么flag存在：
$ samtools view 293T_k18ac_rep1.sam|awk -F '\t' '{print $2}'|sort|uniq -c >PE_flag.statistics
$ cat PE_flag.statistics
43064 113
2144 117
6849 121
45060 129
6769 133
3410 137
117398 141
271561 145
21818023 147
275516 161
21813785 163
43064 177
6849 181
2144 185
11512 2113
69423 2115
173 2121
27242 2129
70403 2131
26102 2145
70479 2147
11239 2161
70214 2163
210 2169
11499 2177
71591 2179
93 2185
27039 2193
72600 2195
27597 2209
71758 2211
11739 2225
71575 2227
94 2233
45060 65
3410 69
6769 73
117398 77
275516 81
21813785 83
271561 97
21818023 99
一下子flag的种类就变多了！双端测序中比较常见的比对成功的FLAG一般有147，163，83，99
这么多种flag，我们就没法一个个单独验证了，接下来我们就先判断这些flag都包含哪些单独的flag，即对这些flag进行拆分。
网友评论