图片 14

在上一篇的文章里我详细介绍了BAM(SAM/CRAM)的格式和一些需要注意的细节,还说了该如何使用samtools在命令行中对其进行操作。但是很多时候这些操作是不能满足我们的实际需要的,比如统计比对率、计算在某个比对质量值之上的read有多少,或者计算PE比对的插入片段长度分布,甚至需要你根据实际情况编写一个新的变异检测算法等。这个时候往往难以直接通过samtools来实现【注】,而是需要编写专门的程序进行计算。因此,在这一篇文章里我们就一起来学习应该如何在程序中借助Pysam来处理BAM文件。

图片 1

【注】关于统计比对率其实是可以通过samtools
stats计算获得的。不过我们这篇文章不是为了争辩samtools能做什么,不能做什么,而是要跟大家讨论该如何编写程序处理BAM。

经过了第四节的长文,我想大家基本上已经知道了一个WGS流程该如何构建起来了吧。但在那一节中限于篇幅有两个很重要的文件我没能展开来讲,分别是:BAM和VCF文件。这篇我们先说BAM文件。

不过,在开始之前我想稍微再补充一下上一节中提到的CRAM——我习惯将其称为BAM的高压缩格式,因为它和BAM/SAM的格式基本相同,但有四点我们需要注意一下:

什么是BAM

BAM是目前基因数据分析中最通用的比对数据存储格式,它既适合于短read也适合于长read,最长可以支持128Mbp的超大read!除了后缀是.bam之外,有些同学可能还会看到.cram,甚至.sam后缀的文件,其实它们一个是BAM的高压缩格式(.cram)——IO效率比原来的BAM要略差;另一个是BAM的纯文本格式(.sam)。当然格式都是一样的,因此为了描述上的清晰,我下面都统一用BAM。

  • CRAM的高压缩是通过借助参考序列和对其他信息的进一步编码来实现的,它相比于BAM有着更高的压缩率,能够节省30%-50%的空间;
  • CRAM目前的IO效率没有BAM高(压得密嘛),约慢30%,但在不断进步,现在已经更新到了3.x版本了;
  • CRAM和BAM可以通过samtools或者picard方便地实现互转;
  • CRAM一定会取代BAM,这话并不是我说的,而是baw/samtools的作者lh3说的。

BAM文件格式

其实一开始它的名字是SAM(The Sequencing Alignment/Map
Format的英文简称),第一次出现的时候,它是bwa比对软件的标准输出文件,但原生的SAM是纯文本文件,十分巨大(比如:一个人30x全基因组测序的sam大小超过600G),非常容易导致存储空间爆满!为了解决这个问题,bwa的开发者李恒(lh3)设计了一种比gz更加高效的压缩算法,对其进行压缩,这就是我们所说的BAM,它的文件大小差不多只有原来的1/6。

在2007年,NGS技术刚刚兴起之时,各类短序列比对软件层出不穷,输出格式也是各有特点,各家各有一套,并没有什么真正的标准可言,可以说那是一个谁都说我最好的时期。

图片 2

但逐渐的,研究者们发现BAM格式对Mapping信息的记录是最全面的,用起来也是最灵活的。bwa的作者还为BAM文件开发了一个非常好用的工具包——Samtools,使得人们对BAM文件的处理变得十分便利,拓展性也变得非常强,后来还有类似于IGV等专门支持BAM的工具也越来越多,因此它就逐渐成为了主流。

现在基本上所有的比对数据都是用BAM格式存储的,俨然已经成为了业内的默认标准。

在2013年,研究者们还专门将Samtools的处理核心剥离出来,并将其打包成为一个专门用于处理高通量数据的API——htslib,除了C语言版本之外还有Java和Python版本,这些在github上都能直接找到。后续许多与NGS数据处理有关的工具基本都会使用这个API进行相关功能的开发,可见其影响力。

ok,背景的介绍就先到此为止了,我们回归主题。下面这个图是我从一份刚刚完成比对的bam文件中截取出来的内容:

图片 3

BAM file

由于屏幕所限,无法把全部的内容都包含进来,特别是header信息,贴在这里仅是为了让还没见过BAM文件的同学们能够对它有一个总体的感觉。

如果是SAM文件,同时你也熟悉linux操作的话,直接在linux终端用less打开即可(注意:不要试图在本地使用文本编辑器,如vim等直接打开文件,会撑死机子的),但如果我们要查看的是BAM,那么必须通过Samtools(可以到samtools的网站下载并安装)。

$ less -SN in.sam          # 打开sam文件
$ samtools view -h in.bam  # 打开bam文件

BAM文件分为两个部分:header和record。这里额外说一句,许多NGS组学数据的存储格式都是由header和record两部分组成的。

以上例子,在samtools view中加上-h参数目的是为了同时把它的header输出出来,如果没有这个参数,那么header默认是不显示的。

图片 4

BAM Header

header内容不多,也不会太复杂,每一行都用‘@’
符号开头,里面主要包含了版本信息,序列比对的参考序列信息,如果是标准工具(bwa,bowtie,picard)生成的BAM,一般还会包含生成该份文件的参数信息(如上图),@PG标签开头的那些。这里需要重点提一下的是header中的@RG也就是Read
group信息,这是在做后续数据分析时专门用于区分不同样本的重要信息。它的重要性还体现在,如果原来样本的测序深度比较深,一般会按照不同的lane分开比对,最后再合并在一起,那么这个时候你会在这个BAM文件中看到有多个RG,里面记录了不同的lane,甚至测序文库的信息,唯一不变的一定是SM的sample信息,这样合并后才能正确处理

其实,关于这一点我在上一篇文章中(引用第四节)讲序列比对时的也特意强调了这些方面,不记得的同学们也可以翻看上一篇的相关内容。

接下来重点要说的是BAM的核心:record(有时候也叫alignment
section,即,比对信息)。这是我们通常所说的序列比对内容,每一行都是一条read比对信息,它的记录看起来是这样的:

图片 5

BAM每一行的内容格式

我这里借用了网上的一张图片来辅助说明,recoed中的每一个信息都是用制表符tab分开的。

下面我们就来仔细瞧瞧这里的每一个信息分别都是什么。

图片 6

BAM格式

以上,前11列是所有BAM文件中都必须要有的信息,而且从描述中我们也能够比较清楚地知道其所代表的含义。但其中,有几个信息实在太重要了,以至于我认为有必要对其进行详细说明。

什么是Pysam

Pysam是一个专门用来处理(BAM/CRAM/SAM)比对数据和变异数据(VCF和BCF)的Python包。它的核心是htslib——一个高通量数据处理API(来自samtools和bwa的核心,基于C语言),开发者们用Python对它直接进行轻量级包装,因此能够在Python中方便地进行调用,并且保证了它与原生C-API功能上的高度一致。

第一,Flag信息

这是一个非常特别并且重要的数字,也是一个容易被忽视的数字,这可能和许多生信工程师也并不完全理解这个值有关。许多同学在第一次看到其官方文档中的描述之后依然会觉得十分困惑,但它里面实际上记录了许多有关read比对情况的信息。想要读懂它的一个关键点是我们不能够将其视为一个数字,而是必须将其转换为一串由0和1组成的二进制码,这一串二进制数中的每一个位(注意是“位”,bit的意思)都代表了一个特定信息,它一共有12位(以前只有8位),所以一般会用一个16位的整数来代表,这个整数的值就是12个0和1的组合计算得来的,因此它的数值范围是0~2048(2的12次方,计算机科学的同学对这种计算应该不陌生)。

那么下面我就结合其文档和自己的实践经验对这12个位的含义用更加通俗易懂的语言来重新描述,如下表:

图片 7

FLAG的含义

为什么是Pysam

因为Pysam可以说是最为官方的版本,有比较固定的开发者在维护,它的稳定性和可靠性都很高。虽然还有一些其它的包同样能够处理BAM但其实它们大多绕不开对htslib的使用,但却没有pysam周全。而且Pysam还集成了tabix的接口,所以除了比对数据之外,还能够用于处理所有用tabix构建过索引的文件,总之就是全且可靠。

如果是文本格式的sam的话,其实也可以直接将其当作普通文本文件来处理,不需借助任何程序包(这在早期的数据分析中经常看到这种操作),只是要麻烦很多(必须自己在程序中处理所有细节,包括解析FLAG和CIGAR信息,以前我也干过不少类似的事情),甚至我还看到有人直接在程序中调用samtools
view把BAM转换成SAM之后再处理的。。。这样的做法实在不推荐。

所以,只要你用的是Python,那么Pysam真的是目前看来比较好的选择。当然如果你用C/C++那么直接用htslib或者bamtools,如果是Java,那么直接使用htsjdk——htslib的java版本。

所以,通过上面这个表的信息,我们就可以清楚地知道每一个FLAG中都包含了什么信息。比如看到FLAG

77时,我们第一步要做的就是将其分解为二进制序列(也可以理解为分解成若干个2的n次方之和):

77 = 000001001101 = 1 + 4 + 8 +64,这样就得到了这个FLAG包含的意思:PE
read,read比对不上参考序列,它的配对read也同样比不上参考序列,它是read1。

当然,如果你希望自己在程序中写一段处理FLAG的代码,那么显然是不会像我们这个例子那样去分解这个整数的,多麻烦啊!那么该如何做呢?其实也很简单,比如我们
要获得其中某个位(假设第N位)的值——只需要将这个FLAG值和2的N次方做与的运算即可。在与运算时,FLAG值首先会被转换成一串二进制序列(如77=000001001101),而2的N次方除了第N位是1之外,其它的都是0,“与”了之后其它信息就会被屏蔽掉。比如,我们想知道该read是否比对上了参考序列,那么只需要计算FLAG
& 4 的值就行了,如果结果是1那么就是比对上了,如果是0则代表没有比上。

不过,在实际工作中,除非遇到特殊的情况,否则我一般更推荐调用官方的htslib这个包来协助处理,它是一个C语言库,如果你用Python,则是pysam——htslib的python包(Java则是htsjdk),包中已经帮我们做了这些处理,可以直接得到结果,下一篇文章里我会用pysam举例说明如何用它来操作bam文件。

另外,下面这一段代码是htslib(samtools的核心库)中定义的12个与flag值进行与操作获取对应位信息的变量,感兴趣的同学可以再htslib里面的sam.h文件中找到,在做一些需要触达基础性原理的开发时或许你会用到。

 /*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
#define BAM_FPAIRED        1
/*! @abstract the read is mapped in a proper pair */
#define BAM_FPROPER_PAIR   2
/*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
#define BAM_FUNMAP         4
/*! @abstract the mate is unmapped */
#define BAM_FMUNMAP        8
/*! @abstract the read is mapped to the reverse strand */
#define BAM_FREVERSE      16
/*! @abstract the mate is mapped to the reverse strand */
#define BAM_FMREVERSE     32
/*! @abstract this is read1 */
#define BAM_FREAD1        64
/*! @abstract this is read2 */
#define BAM_FREAD2       128
/*! @abstract not primary alignment */
#define BAM_FSECONDARY   256
/*! @abstract QC failure */
#define BAM_FQCFAIL      512
/*! @abstract optical or PCR duplicate */
#define BAM_FDUP        1024
/*! @abstract supplementary alignment */
#define BAM_FSUPPLEMENTARY 2048

如何使用Pysam

Pysam

首先,要为我们的Python环境安装这个包,如果已安装过的话可以忽略这一步。有两个方法,pip和bioconda,都比较简单,我们这里以pip——Python的包管理工具来进行:

$ pip install pysam

安装完成之后我们就可以在Python程序中调用pysam了。

第二,CIGAR

CIGAR是Compact Idiosyncratic Gapped Alignment
Report的首字母缩写,称为“雪茄”字符串。

图片 8

作为一个字符串,它用数字和几个字符的组合形象记录了read比对到参考序列上的细节情况,读起来要比FLAG直观友好许多,只是记录的是不同的信息。比如,一条150bp长的read比对到基因组之后,假如看到它的CIGAR字符串为:33S117M,其意思是说在比对的时候这条read开头的33bp在被跳过了(S),紧接其后的117bp则比对上了参考序列(M)。这里的S代表软跳过(Soft
clip),M代表匹配(Match)。CIGAR的标记字符有“MIDNSHP=XB”这10个,分别代表read比对时的不同情况:

图片 9

CIGAR的含义

除了最后‘=XB’非常少见之外,其它的标记符通常都会在实际的BAM文件中碰到。另外,对于M还是再强调一次,CIGAR中的M,不能觉得它代表的是匹配就以为是百分百没有任何miss-match,这是不对的,多态性碱基或者单碱基错配也是用M标记!

读取BAM/CRAM/SAM文件

Pysam中的函数有很多,但是重要的读取函数主要有:

  • AlignmentFile:读取BAM/CRAM/SAM文件
  • VariantFile:读取变异数据(VCF或者BCF)
  • TabixFile:读取由tabix索引的文件;
  • FastaFile:读取fasta序列文件;
  • FastqFile:读取fastq测序序列文件;

等以上几个,其中尤以AlignmentFile和VariantFile为核心。需要我们注意到的地方是,Pysam中的有些函数由于历史原因存在重复,比如名字上只有大小写的差异,但功能却是一样的(比如下图的TabixFile),有些则只是简化了函数名,这些情况用的时候留个心眼就行了。

Tabixfile重名

另外,这篇文章的目的是介绍如何处理比对文件,所以我打算只介绍AlignmentFile

读取比对文件前,我建议先使用samtools index为比对文件构建好索引。当然如果是SAM文件就不必了——它是文本文件,索引的作用是让我们可以对文件进行随机读取,而不必总是从头开始。

下面我先用一个例子作为引子:

import pysam
bf = pysam.AlignmentFile('in.bam', 'rb')

我在这个例子里面,先在程序中导入pysam包,然后调用AlignmentFile函数读取in.bam文件,并把句柄赋值给了bf,bf其实是一个迭代器——Python中的术语,意思就是适合在for循环中进行遍历的对象。

这样我们就是可以通过bf获取这份比对文件中的内容了。比如我们想把in.bam中每一条read的比对位置(包含染色体编号和位置信息),比对质量值和插入片段长度输出(我们的in.bam来自PE测序数据的结果),那么可以这样做:

import pysam
bf = pysam.AlignmentFile('in.bam', 'rb')
for r in bf:
  print r.reference_name, r.pos, r.mapq, r.isize

是不是很简单!打开in.bam文件之后,用for循环对其从头到尾地遍历,并把每个值都赋给r,r在这里代表的就是比对的read信息,它是一个对象(在Pysam由AlignedSegment定义),通过它就可以获取所有的比对信息,比如上面例子中:

  • r.reference_name代表read比对到的参考序列染色体id;
  • r.pos代表read比对的位置;
  • r.mapq代表read的比对质量值;
  • r.isize代表PE read直接的插入片段长度,有时也称Fragment长度;

这里例子的结果如下:

chrM 160 50 235
chrM 161 30 -283
chrM 314 60 -207
...

另外,由于bf是一个迭代器,我们其实还可以用bf.next()一个一个地对其进行访问,而不必在for循环中遍历,这在一些特殊的情况下,这个做法是非常有用且方便的。

当然,上面这个例子其实非常简单,实际上r变量中还有很多其它关于比对的信息,下面这个截图,就是变量中能够获取到的所有比对相关的信息,有好几十个。

函数很多

眼尖的同学可能也发现了,这里面存在一些名字类似的变量,如:r.mapping_quality

r.mapq,它们其实都是比对质量值。类似的也还有几个,这都是上面我提到的历史原因所致,不过这种多余变量随着Pysam的维护也正在逐步变少。

此外,Pysam中的位点坐标体系是0-base(意思是染色体的起始位置是从0而不是1开始算的)而不是1-base,所以上面的输出的160,其实真实位置应该要+1,也就是161

还有,上文我也说过,AlignmentFile除了能够读/写BAM之外,还同样能够读/写CRAM和SAM。区别就在于函数中的第二个参数,比如上面例子中的字符’b’就是用于明确指定BAM文件,’r’字符代表“只读”模式(read首字母)。如果要打开CRAM文件,只需要把b换成c(代表CRAM)就行了,如下:

import pysam
cf = pysam.AlignmentFile('in.cram', 'rc')

那么,如果是SAM文件呢?去掉b或c即可:

import pysam
sf = pysam.AlignmentFile('in.sam', 'r')

第三,MAPQ,比对质量值

这个值同样非常重要,它告诉我们的是这个read比对到参考序列上这个位置的可靠程度,用错误比对到该位置的概率值(转化为Phred
scale)来描述:$$-10logP{错比概率}$$。

因此MAPQ(mapping
quality)值大于30就意味着错比概率低于0.001(千分之一),这个值也是我们衡量read比对质量的一个重要因子。

剩下的几列在上面的格式表中描述的也比较清楚,基本没有过于隐藏的信息,因此我就不打算再一一细说了,如果大家依然有困惑可以到后台留言。

此外,细心的同学可能也已经发现了:fastq的所有信息都被涵盖到了BAM文件中了,包括比对不上的read也在,因此获得了BAM其实也等于获得了所有的read。而且,fastq有时也会被转换成一种uBam文件,指的就是un-mapping
BAM——没有做过比对的BAM文件。它相比于Fastq可以用metadata存储更多有用的信息,不过这不是我们这篇文章想说的内容。

最后,还是再说明一次:BAM文件中除了必须的前11列信息之外,不同的BAM文件中后面记录metadata的列是不固定的,在不同的处理软件中输出时也会有所不同,我们也可以依据实际的情况增删不同的metadata信息

读取特定比对区域内的数据

有时候我们并不需要遍历整一份BAM文件,我们可能只想获得区中的某一个区域(比如chrM中301-310中的信息),那么这个时候可以用Alignmen模块中的fetch函数:

import pysam
bf = AlignmentFile('in.bam', 'rb')
for r in bf.fetch('chrM', 300, 310):
    print r
bf.close()

通过fetch函数就可以定位特定区域了,非常方便。不过,这个时候输入文件in.bam就必须要有索引,不然无法实现这种读取操作。最后用完了,要记得关闭文件——bf.close()。

使用samtools view查看BAM文件

BAM文件由于是特殊的二进制格式,因此没办法通过文本的形式直接打开,要用samtools的view功能在终端上进行查看(上文也已经说到这里在进行系统补充),如:

$ samtools view in.bam

如果不想从头开始看,希望快速地跳转到基因组的其它位置上,比如chr22染色体,那么可以先用samtools
index生成BAM文件的索引(如果已经有索引文件则不需该步骤),然后这样操作:

$ samtools index in.bam  # 生成in.bam的索引文件in.bam.bai
$ samtools view in.bam chr22            # 跳转到chr22染色体
$ samtools view in.bam chr22:16050103   # 跳转到chr22:16050103位置
$ samtools view in.bam chr22:16050103-16050103  # 只查看该位置

来个稍微难一点的例子

问题:如何输出覆盖在某个位置上,比对质量值大于30的所有碱基?

这个问题包含两个部分:

  • 固定的某个位置(我们这里还是用chrM 301这个位置)
  • read比对质量值必须是大于30

如何做呢?这个时候我们要用AlignmentFile模块的另一个函数——pileups来协助解决,代码如下:

import pysam
bf = pysam.AlignmentFile("in.bam", "rb" )
for pileupcolumn in bf.pileup("chrM", 300, 301):
    for read in [al for al in pileupcolumn.pileups if al.alignment.mapq>30]:
        if not read.is_del and not read.is_refskip:
            if read.alignment.pos + 1 == 301:
                print read.alignment.reference_name,\
                      read.alignment.pos + 1,\
                    read.alignment.query_sequence[read.query_position]

bf.close()

这段代码看起来虽然简单,但其实包含了很多信息。总的来说,就是通过pileup获取了所有覆盖到该位置的read,并将其存到pileupcolumn中。然后,对pileupcolumn调用pileups(注意多了一个s)获得一条read中每个比对位置的信息(一条read那么长,并非只覆盖了一个位置),然后通过判断语句留下覆盖到目标位点(301)的碱基。代码中的read.alignment是Pysam中AlignedSegment对象,它包含的内容和上述其它例子中的r是一样的。read.alignment.pos

  • 1还是0-base的原因。最后结果如下:

    chrM 301 A
    chrM 301 A
    chrM 301 A
    chrM 301 C
    chrM 301 C
    chrM 301 C
    chrM 301 C
    chrM 301 C
    chrM 301 C
    chrM 301 C

IGV或者samtools tview查看比对情况

以上,我基本上列举了我们会在终端上如何查看BAM文件的几个最常用操作。但如果你想更直观查看的BAM文件,IGV是目前最好的一个选择,但仅适合于文件还比较小的情况,效果如下:

图片 10

IGV

如果你的BAM文件很大,都超过了你的本地电脑磁盘了,你还是想看该怎么办?你有两个选择:

第一,把你想查看的那部分区域用samtools
view提取出来,生成一份小一些的BAM,然后下载下来,在导入到IGV中。

$ samtools view -h in.bam chr22:16050103-16050203 | samtools view -Sb - > small.bam 

第二,不下载,直接在终端用samtools tview进行查看。samtools
tview有类似于IGV的功能,虽然体验会稍差一些。

$ samtools tview --reference hg38.fa in.bam  

图片 11

samtools tview

在该模式下,按下键盘‘g’后,会跳出一个Goto框,在里面输入想要调整过去的位置,就行了,比如:

图片 12

tview goto

按下esc键则可以取消。另外,为了节省空间,加快查询效率,read中与参考序列相同的部分被用一串串不同颜色的点表示,只留下miss-match的碱基和发生indel变异的区域。其中圆点表示正链比对,逗号表示负链比对。不同的颜色代表不同的比对质量值:白色>=30,黄色20-29,绿色10-19,蓝色0-9。如果你还想知道的其他的功能,可以在tview模式里按下“?”问号,就会弹出类似下面这样的帮助窗口,然后按照指引做就行了。

图片 13

tiview help

虽然看起来不如IGV体验那样好,功能也比较单一(仅可以查看比对情况),但可贵之处在于可以在终端里面直接操作,当需要快速查看某个位置的比对情况时,操作效率非常高。而如果要退出该模式,也非常简单,按下q键就可以了。

创建BAM/CRAM/SAM文件

最后这个例子,我想告诉大家该如何用Pysam输出BAM/CRAM/SAM格式,具体还是看代码吧,这里想输出结果是BAM文件,所以输出模式是“wb”,例子中我们只输出一条比对结果作为说明。

import pysam

header = {'HD': {'VN': '1.0'},
          'SQ': [{'LN': 1575, 'SN': 'chr1'},
                 {'LN': 1584, 'SN': 'chr2'}]
}

tmpfilename = "out.bam"
with pysam.AlignmentFile(tmpfilename, "wb", header=header) as outf:
    a = pysam.AlignedSegment()  # 定义一个AlignedSegment对象用于存储比对信息
    a.query_name = "read_28833_29006_6945"
    a.query_sequence="AGCTTAGCTAGCTACCTATATCTTGGTCTTGGCCG"
    a.flag = 99
    a.reference_id = 0
    a.reference_start = 32
    a.mapping_quality = 20
    a.cigar = ((0,10), (2,1), (0,25))
    a.next_reference_id = 0
    a.next_reference_start=199
    a.template_length=167
    a.query_qualities = pysam.qualitystring_to_array("<<<<<<<<<<<<<<<<<<<<<:<9/,&,22;;<<<")
    a.tags = (("NM", 1),
              ("RG", "L1"))

    outf.write(a)

小结

那么,有关BAM格式的内容我们就暂且先到这里吧,大家如果有疑惑或者感兴趣的内容都可以到后台留言,我都会定时进行回复。在下一篇文章中,我们将重点介绍如何使用pysam来操作bam文件了。


本文首发于我的个人公众号:解螺旋的矿工,欢迎关注,更及时了解更多信息

图片 14

解螺旋的矿工

小结

我写这篇文章的目的主要有两个:第一,充实上一篇文章中关于如何操作BAM的内容;第二,介绍Pysam这一个值得使用的包给大家。另外,我上面列举的例子其实都比较偏于基础操作,这可能和我自身对认知的看法有关。我一直认为,只有真正理解并灵活地应用基础操作,才可以灵活地解决一切复杂的问题。

而且,上面几个例子中用到的模块和函数其实都是比较常用的,所以我比较推荐优先掌握它们。这些例子里面用到的数据我已放到github了,感兴趣的同学可以在公众号后台回复“WGS”即可获得,后续也会陆续有其它的代码和数据可供参考。

最后,Pysam的内容其实还有很多,我所介绍的也仅在对于比对数据的处理,其它很多的模块和函数,包括对Fasta,Fastq,VCF,BCF和Tabix文件的处理,我就不进行一一介绍了,建议大家在使用的时候多看看它的完整文档


本文首发于我的个人公众号:解螺旋的矿工,欢迎关注,更及时了解更多信息

解螺旋的矿工

admin

相关文章

发表评论

电子邮件地址不会被公开。 必填项已用*标注