查看: 2590|回复: 5

[流程分析] GATK4.0和全基因组数据分析实践(上)

[复制链接]
  • TA的每日心情
    害羞
    2019.9.24 10:09
  • 签到天数: 33 天

    连续签到: 1 天

    [LV.5]常住居民I

    迅猛龙

    Rank: 8Rank: 8

    主题
    161
    奥币
    1940
    积分
    1899
    注册时间
    2015.11.14
    在线时间
    78 小时

    发表于 2019.8.12 11:11:37 | 显示全部楼层 |阅读模式
    本帖最后由 犬夜× 于 2019.8.12 11:18 编辑

    写在前面:本文转载自微信公众号:碱基矿工,是WGS系列文章中的第4节 构建WGS主流程中的扩展和延伸。感兴趣的同学可以了解之前的整个WGS系列,保证让你收获满满。

    下面是正文:
    ======================================

    本次文章分为上下两篇,这是第一篇,都是WGS系列第四节的延伸和补充,本意是结合具体的数据,让更多的人能够更好地理解整个WGS数据的分析和处理过程,我也结合自身的工作经验给出一些做项目过程中的建议,以作参考,希望能够对你有更多的帮助。另外,接下来我将系统写一个关于全基因组关联分析(GWAS)的文章,同时还会有更多全面而且紧扣前沿的技术文章分享出来。

    那么,事不宜迟我们马上开始。考虑到实际的数据分析环境,我们只在Linux命令行终端(Terminal)中进行,执行步骤都写在shell,因此不会有窗口式的操作(使用Mac OS的同学可以使用Mac自带的Terminal,与Linux操作一致),在这篇文章中我们会用到以下几个工具:
    • sratoolkit

    • bgzip

    • tabix

    • bwa

    • samtools

    • GATK 4.0


    这些软件都可以在github上找到(包括GATK),需要各位自行安装。这里补充一句,目前GATK4.0的正式版本已经发布,它的使用方式与之前相比有着一些差异(变得更加简单,功能也更加丰富了),增加了结构性变异检测和很多Spark、Cloud-Only的功能,并集成了MuTect2和picard的所有功能(以及其他很多有用的工具),这为我们减少了许多额外的工具,更加有利于流程的构建和维护,4.0之后的GATK是一个新的篇章,大家最好是掌握这一个版本!另外,3.x的版本貌似也已经不提供下载通道了,如果你还想使用3.x的话可以在公众号后台回复“GATK3”,我为你准备了一个GATK官方3.7的版本。我们这里则使用最新的4.0版本。

    项目目录结构
    清晰的目录结构是管理众多项目的有效途径,经久不忘,随时可查。虽然看起来有些原始,但在Linux终端下面,我目前还没有发现更好的文件管理办法。这个项目的目录结构,我的建议是按照时间+项目的规则来命名,下面是我的目录结构:
    [AppleScript] 纯文本查看 复制代码
    ./201802_wgs_practice/
    ├── bin
    ├── input
    └── output
    顶层的项目名就是20180203_wgs_practice,下面有三个主目录:
    • input:存储所有输入数据

    • output:存储所有输出数据

    • bin:存放所有执行程序和代码


    output只存放结果数据,它是由input和bin中的数据和程序流程生成的。这样做的好处是层次分明,流程逻辑清楚,数据互不干扰。

    使用E.coli K12完成比对和变异检测

    人类基因组数据很大,参考序列长度是3Gb。而一个人的高深度测序数据往往是这个数字的30倍——100Gb。如果直接用这样的数据来完成本文的分析,那么许多同学需要下载大量的原始数据。除了下载时间很长之外,如果没有合适的集群,只是在自己的桌面电脑上干这样的事情,那么硬盘空间也将很快不够用。而且,要在单机电脑上完成这样一个高深度WGS数据的分析,除了对机器性能有要求之后,跑起来也需要花上差不多140个小时——相信大家都等不起呀。

    因此,为了解决这个问题,我找了E.coli K12(一种实验用的大肠杆菌)的数据作为代替,用来演示数据比对变异检测两个最消耗计算资源和存储空间的步骤。E.coli K12的特点是数据很小,它的基因组长度只有4.6Mb,很适合大家用来快速学习WGS的数据分析,遇到人类的数据时,再做替换就行了。

    下载E.coli K12的参考基因组序列

    熟悉的同学应该第一时间能够知道,这些物种的基因组参考序列都可以在NCBI上获取,我们这里也是一样,可以在NCBI网站上直接搜索这个序列,为了简化步骤,我直接给出E.coli K12参考序列的ftp地址给大家下载之用:

    [AppleScript] 纯文本查看 复制代码
    ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz

    你可以在Linux(或者Mac OSX)命令行上直接使用wget,将这个fasta下载下来,由于它很小,所以几秒之后我们就可以得到这个fasta序列。
    [AppleScript] 纯文本查看 复制代码
    $ wget [url=ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz]ftp://ftp.ncbi.nlm.nih.gov/genom ... 84v2_genomic.fna.gz[/url]

    为了接下来表达上的清晰和操作上的方便,我们使用bgzip将这个序列文件进行解压并把名字重命名为E.coli_K12_MG1655.fa,这样就一目了然了。
    [AppleScript] 纯文本查看 复制代码
    $ gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa

    E.coli K12只有一条完整的染色体,你打开文件后将会看到和我一样的内容:
    [AppleScript] 纯文本查看 复制代码
    >NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
    AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
    GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC
    接着,我们用samtools为它创建一个索引,这是为方便其他数据分析工具(比如GATK)能够快速地获取fasta上的任何序列做准备。
    [AppleScript] 纯文本查看 复制代码
    $ /Tools/common/bin/samtools faidx E.coli_K12_MG1655.fa


    这时会生成一份E.coli_K12_MG1655.fa.fai文件。除了方便其他工具之外,我们可以通过这样的索引来获取fasta文件中任意位置的序列或者任意完整的染色体序列。可以很方便地完成对参考序列(或者任意fasta文件)特定区域序列的提取。举个例子:

    [AppleScript] 纯文本查看 复制代码
    $ samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200
    我们就获得了E.coli K12参考序列上的这一段序列:
    [AppleScript] 纯文本查看 复制代码
    >NC_000913.3:1000000-1000200
    GTGTCAGCTTTCGTGGTGTGCAGCTGGCGTCAGATGACAACATGCTGCCAGACAGCCTGA
    AAGGGTTTGCGCCTGTGGTGCGTGGTATCGCCAAAAGCAATGCCCAGATAACGATTAAGC
    AAAATGGTTACACCATTTACCAAACTTATGTATCGCCTGGTGCTTTTGAAATTAGTGATC
    TCTATTCCACGTCGTCGAGCG
    这个小技巧在特定的时候非常实用

    下载E.coli K12的测序数据

    基因组参考序列准备好之后,接下来我们需要下载它的测序数据。E.coli K12作为一种供研究使用的模式生物,自然已经有许多的测序数据在NCBI上了,在这里我们选择了其中的1个数据——SRR1770413。这个数据来自Illumina MiSeq测序平台(不用担心平台的事情),read长度是300bp,测序类型Pair-End(没了解过PE read同学可以参考我前面WGS系列第四节文章)。你可以在NCBI上直接搜到:https://www.ncbi.nlm.nih.gov/sra/?term=SRR1770413

    在NCBI给出的信息页面中,我们可以清楚地看到这个数据的大小(如下图)——差不多200MB,一般家庭网速也能够较快下载完成。


    从NCBI上下载下来的测序数据,不是我们熟悉的fastq格式,而是SRA(一种NCBI自己设计的测序数据存储格式,具较高的压缩率),我们需要对其进行转换,下文详述。现在我们先下载,有两个下载方式(我在这里告诉大家的方法同样适用于其他类型的数据),第一个是如上面所说搜索到SRR1770413这个数据的ftq地址,然后直接在命令行中执行wget进行下载,如下:

    [AppleScript] 纯文本查看 复制代码
    $ wget [url=ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra]ftp://ftp-trace.ncbi.nlm.nih.gov ... 0413/SRR1770413.sra[/url]

    注意,下载下来的这个SRA文件虽然只有一份,但是里面其实存了read1和read2的测序数据,我们要把它解出来,转换成为我们所需的fastq格式。这个时候,我们就需要用到NCBI的官方工具包sratoolkit—— https://github.com/ncbi/sra-tools/wiki/Downloads ,大家下载对应系统的版本,直接解压之后就可以使用了。


    sratoolkit是一个工具包,所有的执行程序都在它解压后的bin文件夹下,我们要把SRA转换为fastq,只需直接通过工具包中的fastq-dump即可完成。
    [AppleScript] 纯文本查看 复制代码
    $ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413.sra
    然后我们就会得到这个E.coli K12数据的read1和read2了:
    SRR1770413_1.fastq
    SRR1770413_2.fastq
    另一个数据下载的方法就是用上面说到的sratoolkit,我也比较推荐这个方法,操作很简单。同样是使用fastq-dump(没错,与上面一样,区别在于下载的时候,输入的是数据的SRA编号),它可以在我们下载的过程中就直接将SRA转换为两个fastq。
    [AppleScript] 纯文本查看 复制代码
    $ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413
    下载完成后,我们最好用bgzip(不推荐gzip)将其压缩为.gz文件,这样可以节省空间,而且也不会对接下来的数据分析产生影响。
    [AppleScript] 纯文本查看 复制代码
    $ /Tools/common/bin/bgzip -f SRR1770413_1.fastq
    $ /Tools/common/bin/bgzip -f SRR1770413_1.fastq

    至此,E.coli K12相关的数据我们就都准备好了,先看一眼我现在的目录结构。

    [AppleScript] 纯文本查看 复制代码
    ./201802_wgs_practice/
    ├── bin
    ├── input
    │   ├── E.coli
    │       ├── fasta
    │       │   ├── E.coli_K12_MG1655.fa
    │       │   ├── E.coli_K12_MG1655.fa.fai
    │       │   └── work.log.sh
    │       └── fastq
    │           ├── SRR1770413_1.fastq.gz
    │           ├── SRR1770413_2.fastq.gz
    │           └── work.log.sh   
    ├── output
    └── work.log.sh

    其中各个目录下的work.log.sh记录了我在该目录下的所有重要操作——这是我的个人习惯,目的是方便以后反查数据的需要。

    数据准备完毕之后,接下来就可以进行具体的分析了。


    本帖子中包含更多资源

    您需要 登录 才可以下载或查看,没有帐号?立即注册

    x
    我的意中人是一个盖世英雄,有一天他会踩着七彩祥云来娶我……
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    2019.9.24 10:09
  • 签到天数: 33 天

    连续签到: 1 天

    [LV.5]常住居民I

    迅猛龙

    Rank: 8Rank: 8

    主题
    161
    奥币
    1940
    积分
    1899
    注册时间
    2015.11.14
    在线时间
    78 小时

     楼主| 发表于 2019.8.12 11:16:56 | 显示全部楼层
    质控

    质控是必须做的,我们需要完整认识原始的测序数据质量到底如何,该步骤不能省略。我专门为此单独写了一篇文章(WGS系列第三节),在正式的数据分析过程中,大家可以参考它来完成数据的质控,然后再进行接下来的分析。本篇文章为了控制篇幅和尽可能扣住核心内容,就不再对此深入展开,大家如果碰到问题,可以到在后台留言,或者加入我的交流圈(知识星球:解螺旋技术交流圈)和更多有经验的人一起交流。

    比对

    首先是比对。所谓比对就是把测序数据定位到参考基因组上,确定每一个read在基因组中的位置。这里,我们依然用目前使用最广的BWA来完成这个工作。在正式比对之前,需要先为参考序列构建BWA比对所需的FM-index(比对索引)。


    [AppleScript] 纯文本查看 复制代码
    $ /Tools/common/bin/bwa index E.coli_K12_MG1655.fa
    
    [bwa_index] Pack FASTA... 0.06 sec
    [bwa_index] Construct BWT for the packed sequence...
    [bwa_index] 1.78 seconds elapse.
    [bwa_index] Update BWT... 0.04 sec
    [bwa_index] Pack forward-only FASTA... 0.03 sec
    [bwa_index] Construct SA from BWT and Occ... 0.63 sec
    [main] Version: 0.7.16a-r1181
    [main] CMD: /Tools/common/bin/bwa index E.coli_K12_MG1655.fa
    [main] Real time: 3.266 sec; CPU: 2.550 sec

    由于这个序列很短,只需几秒就可以完成这个索引文件的构建(对于人类基因组则需要3个小时的时间)。创建完毕之后,将多出5份以E.coli_K12_MG1655.fa为前缀的序列索引文件。
    E.coli_K12_MG1655.fa.amb
    E.coli_K12_MG1655.fa.ann
    E.coli_K12_MG1655.fa.bwt
    E.coli_K12_MG1655.fa.pac
    E.coli_K12_MG1655.fa.sa
    现在我们使用bwa完成比对,用samtools完成BAM格式转换、排序并标记PCR重复序列。步骤分解如下:
    [AppleScript] 纯文本查看 复制代码
    #1 比对
    time /Tools/common/bin/bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_1.fastq.gz /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_2.fastq.gz | /Tools/common/bin/samtools view -Sb - > /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** bwa mapping done **"#2 排序
    time /Tools/common/bin/samtools sort -@ 4 -m 4G -O bam -o /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** BAM sort done"rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
    
    #3 标记PCR重复
    time /Tools/common/bin/gatk/4.0.1.2/gatk MarkDuplicates -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -M /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup_metrics.txt && echo "** markdup done **"#4 删除不必要文件(可选)
    rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
    rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam
    
    #5 创建比对索引文件

    从上面的命令大家也可以看到,我严格按照上文提到的项目目录规范来执行(步骤中涉及到的数据路径也都尽可能使用全路径),这个比对的shell存放在bin目录下,名称是bwa_and_markdup.sh,从名称也能够一眼可以看出是要做什么的。

    简单解释一下这个shell所做的事情:首先利用bwa mem比对模块将E.coli K12质控后的测序数据定位到其参考基因组上(我们这里设置了4个线程来完成比对,根据电脑性能可以适当调大),同时通过管道('|' 操作符)将比对数据流引到samtools转换为BAM格式(SAM的二进制压缩格式),然后重定向('>'操作符)输出到文件中保存下来。

    -R 设置Read Group信息,虽然我在以前的文章中已经反复强调过它的重要性,但这里还是再说一次,它是read数据的组别标识,并且其中的ID,PL和SM信息在正式的项目中是不能缺少的(如果样本包含多个测序文库的话,LB信息也不要省略),另外由于考虑到与GATK的兼容关系,PL(测序平台)信息不能随意指定,必须是:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这12个中的一个。

    接着用samtools对原始的比对结果按照参考序列位置从小到大进行排序(同样是4个线程),只有这个步骤完成之后才可以继续往下。

    然后,我们使用GATK标记出排完序的数据中的PCR重复序列。这个步骤完成后,如无特殊需要,我们就可以直接删除前面那两个BAM文件了(原始比对结果和排序后的结果)——后续几乎不会再用到那两份文件了。关于标记PCR重复序列的操作比较简单,不再细说(如果希望了解更多有关重复序列特征的信息可以回看WGS系列第四节中的内容)。

    最后,我们再用samtools为E_coli_K12.sorted.markdup.bam创建索引。我认为不论是否有后续分析,为BAM文件创建索引应该作为一个常规步骤,它可以让我们快速地访问基因组上任意位置的比对情况,这一点非常有助于我们随时了解数据。

    至于每个步骤最前面的time,则是用于记录执行时间的,有助于我们清楚地知道每一个分析过程都花了多少时间,当需要优化流程的时候这个信息会很有用。

    变异检测

    接下来是用GATK完成变异检测。但在开始之前之前我们还需要先为E.coli K12的参考序列生成一个.dict文件,这可以通过调用CreateSequenceDictonary模块来完成(这是原来picard的功能)。
    [AppleScript] 纯文本查看 复制代码
    $ /Tools/common/bin/gatk/4.0.1.2/gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"
    唯一需要注意的是.dict文件的名字前缀需要和fasta的一样,并跟它在同一个路径下,这样GATK才能够找到。

    OK,现在我们就可以进行变异检测了,同样使用GATK 4.0的HaplotypeCaller模块来完成。由于我们只有一个样本,要完成这个工作其实很简单,直接输入比对文件和参考序列就行了,但是考虑到实际的情况,我想告诉大家一个更好的方式(虽然这会多花些时间),就是:先为每个样本生成一个GVCF,然后再用GenotypeGVCFs对这些GVCF进行joint calling,如下 ,我把命令都写在gatk.sh中,并执行。
    [AppleScript] 纯文本查看 复制代码
    #1 生成中间文件gvcf
    time /Tools/common/bin/gatk/4.0.1.2/gatk HaplotypeCaller \
     -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
     --emit-ref-confidence GVCF \
     -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam \
     -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf && echo "** gvcf done **"
    
    #2 通过gvcf检测变异
    time /Tools/common/bin/gatk/4.0.1.2/gatk GenotypeGVCFs \
     -R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
     -V /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf \
     -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf && echo "** vcf done **"
    很快我们就获得了E.coli K12这个样本初步的变异结果——E_coli_K12.vcf。之所以非要分成两个步骤,是因为我想借此告诉大家,变异检测不是一个样本的事情,有越多的同类样本放在一起joint calling结果将会越准确,而如果样本足够多的话,在低测序深度的情况下也同样可以获得完整并且准确的结果,而这样的分步方式是应对多样本的好方法。

    最后,我们用bgzip对这个VCF进行压缩,并用tabix为它构建索引,方便以后的分析。

    [AppleScript] 纯文本查看 复制代码
    #1 压缩 
    time /Tools/common/bin/bgzip -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf
    
    #2 构建tabix索引
    time /Tools/common/bin/tabix -p vcf /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf.gz
    bgzip压缩完成之后,原来的VCF文件会被自动删除。

    为了保持一致,现在再看一下完成到这里之后我们的目录长什么样了,供大家对照。
    [AppleScript] 纯文本查看 复制代码
    ./201802_wgs_practice/
    ├── bin
    │   ├── bwa_and_markdup.sh
    │   └── gatk.sh
    ├── input
    │   └── E.coli
    │       ├── fasta
    │       │   ├── E.coli_K12_MG1655.dict
    │       │   ├── E.coli_K12_MG1655.fa
    │       │   ├── E.coli_K12_MG1655.fa.amb
    │       │   ├── E.coli_K12_MG1655.fa.ann
    │       │   ├── E.coli_K12_MG1655.fa.bwt
    │       │   ├── E.coli_K12_MG1655.fa.fai
    │       │   ├── E.coli_K12_MG1655.fa.pac
    │       │   ├── E.coli_K12_MG1655.fa.sa
    │       │   └── work.log.sh
    │       └── fastq
    │           ├── SRR1770413_1.fastq.gz
    │           ├── SRR1770413_2.fastq.gz
    │           └── work.log.sh
    ├── output
    │   └── E.coli
    │       ├── E_coli_K12.g.vcf
    │       ├── E_coli_K12.g.vcf.idx
    │       ├── E_coli_K12.sorted.markdup.bam
    │       ├── E_coli_K12.sorted.markdup.bam.bai
    │       ├── E_coli_K12.sorted.markdup_metrics.txt
    │       ├── E_coli_K12.vcf
    │       └── E_coli_K12.vcf.idx
    └── work.log.sh

    如果大家仔细看过WGS系列第四节的话,会发现我这里缺少了两个步骤:重比对和BQSR。没有执行BQSR是因为E.coli K12没有那些必须的known变异集(或者有但我没找到),所以无法进行;但没有重比对,则是因为我在GATK 4.0中没发现IndelRealigner这个功能,虽然我们使用GATK HaplotypeCaller或者Mutect2的话确实可以省略这个步骤,但如果是其他软件来进行变异检测那么该步骤依然十分重要,我目前不太清楚为何GATK 4.0没有将这个功能单独分离出来。

    后面要谈到的就是变异的质控了。很遗憾我们这个E.coli K12的变异结果并不适合通过VQSR来进行过滤,原因上面也提到了一些,它不像人类的基因组数据,有着一套适合用来训练过滤模型的已知变异集(dbSNP,1000G,Hapmap和omini等)。其实这种情况有时候我们在工作中也会碰到,比如有些捕获测序(Panel测序数据,甚至外显子测序)的数据,由于它的区域较小,获得的变异也不多,导致最终没法满足VQSR进行模型训练时所需的最低变异数要求,那时你也不能通过这个方式协助变异质控。那么碰到这种情况的时候该怎么办?我将这部分的内容放在了下一篇文章中,在那里我们再来讨论这个问题。我也会告诉大家变异质控的基本逻辑,而不是简单罗列一个命令,同时也会再用NA12878这个人的数据来进一步告诉大家如何比较和评估变异结果。

    小结

    至此,这个篇文章的上半部分就到此为止了。除了那些重要的内容之外,在上文中,你会看到我反复提到了创建“索引”这个事情,比如为fasta,为BAM,为VCF。我为什么非要反复强调这个事情不可呢?因为我发现许多初学者并不知道索引的作用,当被问到如何从巨大的比对文件或者变异文件中提取某个信息时,总是要走弯路——努力写程序去提取,既慢又费力,结果还不一定好,甚至有些有一定经验的同学也不知道使用bgzip和tabix的好处,因此我才反复在文章里提及。



    我的意中人是一个盖世英雄,有一天他会踩着七彩祥云来娶我……
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    害羞
    2020.6.11 15:25
  • 签到天数: 735 天

    连续签到: 1 天

    [LV.9]以坛为家II

    迅猛龙

    Rank: 8Rank: 8

    主题
    26
    奥币
    4338
    积分
    1804
    注册时间
    2016.1.8
    在线时间
    431 小时

    发表于 2019.8.12 13:53:31 | 显示全部楼层
    666
    新的一天加油!
    回复

    使用道具 举报

  • TA的每日心情
    yes!
    2020.5.22 10:21
  • 签到天数: 318 天

    连续签到: 1 天

    [LV.8]以坛为家I

    帝王蝶

    Rank: 4

    主题
    0
    奥币
    549
    积分
    251
    注册时间
    2018.1.22
    在线时间
    116 小时

    发表于 2019.8.14 09:51:00 | 显示全部楼层
    lalalala!
    回复

    使用道具 举报

  • TA的每日心情
    害羞
    2019.9.24 10:09
  • 签到天数: 33 天

    连续签到: 1 天

    [LV.5]常住居民I

    迅猛龙

    Rank: 8Rank: 8

    主题
    161
    奥币
    1940
    积分
    1899
    注册时间
    2015.11.14
    在线时间
    78 小时

     楼主| 发表于 2019.8.15 10:39:07 | 显示全部楼层
    我的意中人是一个盖世英雄,有一天他会踩着七彩祥云来娶我……
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    害羞
    2019.9.24 10:09
  • 签到天数: 33 天

    连续签到: 1 天

    [LV.5]常住居民I

    迅猛龙

    Rank: 8Rank: 8

    主题
    161
    奥币
    1940
    积分
    1899
    注册时间
    2015.11.14
    在线时间
    78 小时

     楼主| 发表于 2019.8.15 10:39:39 | 显示全部楼层
    我的意中人是一个盖世英雄,有一天他会踩着七彩祥云来娶我……
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    快速回复 返回顶部 返回列表