annovar使用笔记

annovar是一款常用的注释软件,可在其官网注册后下载。
annovar无需安装,下载后解压即可直接使用。annovar软件里面是几个perl写的脚本:

annotate.pl 下载数据库,注释数据
retrieve_seq_from_fasta.pl
coding_change.pl 可用来推断蛋白质序列
convert2annovar.pl 将变异文件转化annovar可以使用的文件格式
table_annovar.pl 注释文件,一次可完成多种类型的注释
variants_reduction.pl 可用来更灵活地定制过滤注释流程


1、Download DataBase

annovar提供了许多常用的数据库文件可使用annotate.pl直接下载:

# 这里下载三个数据库 refgene、dbsnp、1000genomes
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar avsnp150 humandb
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar 1000g2015aug humandb

2、使用convert2annovar.pl将输入文件进行格式转换

使用annovar注释对输入文件有一定的格式要求,因此注释前需对输入文件的格式做简单的转换。
annovar对输入文件有明确格式要求的只有前5列,这5列依次必须为: Chromosome ("chr" prefix is optional), Start, End, Reference Allelel, Alternative Allele. 其余的可以列随需要添加。输入文件格式如下:

1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss

其中插入或者删除以-表示, “0” means this information is not readily available.
使用convert2annovar.pl最常用的就是对vcf文件进行转换:

perl convert2annovar.pl -format vcf4 example/ex2.vcf > ex2.avinput

除了这种最简单常用的用法外,convert2annovar.pl还有一些非常有用的参数.

2.1 -allsample

对于含有多个样本的vcf文件,格式转换时只会取其第一个样本进行注释,也就是说即使别的样本在这个位点有变异,只要第一个样本在某个位点没有变异转换时就会将这个位点去掉不会出现在注释文件中。如果想要得到所有样本的变异位点的注释的话,可以先将其拆分为几个样本的注释输入文件:

# 转换格式时vcf中的每一个样本会单独生成一个待注释的vcf文件
perl convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -allsample

这一功能另一用法类似于vcftools中的vcf-subset,可将多样本的vcf文件分开为多个单样本的vcf文件,但是效率要比vcf-subset高得多。

2.2 -includeinfo, -comment

includeinfo参数会保留vcf文件中的所有信息。
comment参数会保留vcf文件的头部注释信息(以#开头的行)。
convert2annovar.pl -format vcf4 example/ex2.vcf -outfile ex2 -allsample -includeinfo -comment

2.3 dbSNP identifiers

有时我们得到了一些有兴趣的snp位点(dbsnp rsID)且只想对这些位点进行注释,可使用-format rsid来完成这一功能:

[kaiwang@biocluster ~/]$ cat example/snplist.txt 
rs74487784
rs41534544
rs4308095
rs12345678

[kaiwang@biocluster ~/]$ convert2annovar.pl -format rsid example/snplist.txt -dbsnpfile humandb/hg19_snp138.txt > snnplist.avinput 
NOTICE: Scanning dbSNP file humandb/hg19_snp138.txt...
NOTICE: input file contains 4 rs identifiers, output file contains information for 4 rs identifiers
WARNING: 1 rs identifiers have multiple records (due to multiple mapping) and they are all written to output

[kaiwang@biocluster ~/]$ cat snplist.avinput 
chr2 186229004 186229004 C T rs4308095
chr7 6026775 6026775 T C rs41534544
chr7 6777183 6777183 G A rs41534544
chr9 3901666 3901666 T C rs12345678
chr22 24325095 24325095 A G rs74487784
2.4 All variants in a genomic region

现在我们发现有一段区域很可疑,其变异可能与我们性状相关,只想对这一段区域进行注释,通样的可以使用-format参数解决:

[kaiwang@biocluster ~/]$ convert2annovar.pl -format region -seqdir humandb/hg19_seq/ chr1:2000001-2000003
NOTICE: Reading region from STDIN ... Done with 1 regions from 1 chromosomes
NOTICE: Finished reading 1 sequences from humandb/hg19_seq/chr1.fa
NOTICE: Finished writting FASTA for 1 genomic regions to stdout
1 2000001 2000001 A C
1 2000001 2000001 A G
1 2000001 2000001 A T
1 2000002 2000002 T A
1 2000002 2000002 T C
1 2000002 2000002 T G
1 2000003 2000003 C A
1 2000003 2000003 C G
1 2000003 2000003 C T

而且,通过一些其他参数还能限定这段区域内的变异类型来进行注释,譬如2bp insertion,3bp的indel等,详见http://annovar.openbioinformatics.org/en/latest/user-guide/input/
此外,annovar还能对MAQ、GFF等格式进行注释。

几个需要注意的地方:
  1. vcf文件在格式转换时,若突变位点有两个不同的等位基因则在结果文件中会分两行放。
  2. 在注释时,遇到格式不符合的行会跳过继续注释而不是终止注释,最后那些格式不符合的行会生成另一个文件(*.invalid_input)。

3、Annotate

annovar提供了两个脚本以供注释使用:annotate_variation.pl一次注释一个数据库,table_annovar.pl一次注释多个数据库。

### 使用annotate_variation.pl注释refgene数据库
perl annotate_variation.pl input.av -buildver hg19 -geneanno -dbtype refGene humandb/  -out ex1 
# -geneanno 表示使用基于基因的注释,另有-filter表示基于过滤的注释,-region表示基于位置的注释
# -dbtype refGene 表示使用"refGene"数据库
# -out ex1 表示输出文件以ex1为前缀,亦可用 --outfile 直接指定文件名

### 使用table_annovar.pl注释多个数据库
perl table_annovar.pl example/ex1.avinput humandb/ -buildver hg19 -out myanno -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f -nastring . -csvout -polish -xref example/gene_xref.txt
#-bulidver hg19 表示使用的参考基因组版本
#-out myanno 表示输出文件前缀,亦可用 --outfile 直接指定文件名
#-remove 表示删除中间文件
#-protocol 表示使用的数据库,其数据库顺序要与后面的operation注释方式对应上
#-operation 表示对应数据库的注释类型(g代表gene-based、r代表region-based、f代表filter-based,gx means gene-based with cross-reference annotation (from -xref argument))
#-nasting . 点号代替缺省值
#-csvout 表示输出为csv格式

4、注释类型

annovar注释类型有三种:gene-based、region-based、filter-based

4.1 gene-based

Gene-based annotation是根据SNPs以及CNVs的位置信息来确定是否会造成编码序列以及开放阅读框的改变从而影响氨基酸的改变,使用者可以自主选择RefSeq genes, 包括UCSC genes, ENSEMBL genes, GENCODE genes, AceView genes等来进行注释。注释后会生成两个文件:ex1.variant_function and ex1.exonic_variant_function。

$ perl annotate_variation.pl -geneanno -dbtype refGene -out ex1 -build hg19 example/ex1.avinput humandb/
$ cat ex1.variant_function
UTR5 ISG15(NM_005101:c.-33T>C) 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
UTR3 ATAD3C(NM_001039211:c.*91G>T) 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
splicing NPHP4(NM_001291593:exon19:c.1279-2T>A,NM_001291594:exon18:c.1282-2T>A,NM_015102:exon22:c.2818-2T>A) 1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
intronic DDR2 1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
intronic DNASE2B 1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
intergenic LOC645354(dist=11566),LOC391003(dist=116902) 1 13211293 13211294 TC - comments: rs59770105, a 2-bp deletion
intergenic UBIAD1(dist=55105),PTCHD2(dist=135699) 1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
intergenic LOC100129138(dist=872538),NONE(dist=NONE) 1 105492231 105492231 A ATAAA comments: rs10552169, a block substitution
exonic IL23R 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
exonic ATG16L1 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
exonic NOD2 16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
exonic NOD2 16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
exonic NOD2 16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
exonic GJB2 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
exonic CRYL1,GJB6 13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss 
第一个文件包括对于所有突变的注释,通过在文件最前面加入两列,以tab分割
第一列为变异所在基因位置的类型,如外显子,内含子,UTR5,UTR3,基因间等
第三列为变异文件原有的comment信息
第二列为对第一列的描述信息,详情见下
detail explain of first column
#ex1.exonic_variant_function
[kaiwang@biocluster ~/]$ cat ex1.exonic_variant_function
line9 nonsynonymous SNV IL23R:NM_144701:exon9:c.G1142A:p.R381Q, 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
line10 nonsynonymous SNV ATG16L1:NM_001190267:exon9:c.A550G:p.T184A,ATG16L1:NM_017974:exon8:c.A841G:p.T281A,ATG16L1:NM_001190266:exon9:c.A646G:p.T216A,ATG16L1:NM_030803:exon9:c.A898G:p.T300A,ATG16L1:NM_198890:exon5:c.A409G:p.T137A, 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
line11 nonsynonymous SNV NOD2:NM_022162:exon4:c.C2104T:p.R702W,NOD2:NM_001293557:exon3:c.C2023T:p.R675W, 16 50745926 50745926 C comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
line12 nonsynonymous SNV NOD2:NM_022162:exon8:c.G2722C:p.G908R,NOD2:NM_001293557:exon7:c.G2641C:p.G881R, 16 50756540 50756540 G comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
line13 frameshift insertion NOD2:NM_022162:exon11:c.3017dupC:p.A1006fs,NOD2:NM_001293557:exon10:c.2936dupC:p.A979fs, 16 50763778 5076377comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
line14 frameshift deletion GJB2:NM_004004:exon2:c.35delG:p.G12fs, 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss line15 frameshift deletion GJB6:NM_001110221:wholegene,GJB6:NM_001110220:wholegene,GJB6:NM_001110219:wholegene,CRYL1:NM_015974:wholegene,GJB6:NM_006783:wholegene, 13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss
第二个输出文件以.exonic_variant_function结尾,只列出外显子(氨基酸会改变)的变异
第一列为第一个文件中该变异所在的行号;
第二列为该变异的功能性后果,如外显子改变导致的氨基酸变化,阅读框移码,无义突变,终止突变等
第三列为基因名称,转录识别标志和相应的转录本的序列变化
第四列为原输入文件内容
4.2 region-based annotation

其与Gene-based annotation作用相反,它是用来确认在特定区域的突变造成的影响。比如在44个物种的保守基因区域,预测的转录因子结合区域,基因重复区域,GWAS分析区域,基因突变数据库,表观组学位点等。此处以Conserved genomic elements annotation为例介绍region-based annotation的使用:

[kaiwang@biocluster ~/]$ annotate_variation.pl -regionanno -build hg19 -out ex1 -dbtype phastConsElements46way example/ex1.avinput humandb/
NOTICE: Reading annotation database humandb/hg19_phastConsElements46way.txt ... Done with 5163775 regions
NOTICE: Finished region-based annotation on 12 genetic variants in ex1.hg19.avinput
NOTICE: Output files were written to ex1.hg19_phastConsElements46way
# -regionanno 表示使用基于区域的注释
# -dbtype phastConsElements46way 表示使用"phastConsElements46way"数据库,注意需要使用Region-based的数据库
######输出文件
[kaiwang@biocluster ~/]$ cat ex1.hg19_phastConsElements46way
phastConsElements46way Score=387;Name=lod=50 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
phastConsElements46way Score=420;Name=lod=68 16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
phastConsElements46way Score=385;Name=lod=49 16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
phastConsElements46way Score=395;Name=lod=54 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
phastConsElements46way Score=545;Name=lod=218 13 20797176 21105944 0 - comments: a 342kb deletion encompassing GJB6, associated with hearing loss
输出文件:输出的注释文件第1列为“phastConsElements46way”,对应注释的类型,这里的phastCons 46-way alignments属于保守的基因组区域的注释;
第二列包含评分和名称,评分来自UCSC,可以使用--score_threshold和--normscore_threshold来过滤评分低的变异,“Name=lod=x”名称表示该区域的名称;
剩余的部分为输入文件的内容
4.3 filter-based annotation

Filter-based annotation是用以确认已记录在特定数据库里的突变。例如想要知道突变是否为novel variation就需要知道该突变是否存在于dbSNP库里,它在1000 genome project里面等位基因频率怎样,以及计算一系列突变项目得分并加以过滤。它区别于region-based annotation就在于它针对突变碱基进行工作,而region-based annotation 针对染色体位置。举例来说就是region-based比对chr1:1000-1000而filter-based比对chr1:1000-1000上的A->G。
它拥有多种数据库,包括针对全基因组测序的突变频率,针对全外显子数据测序的突变频率,在孤立或者小类群人群中的突变频率,全基因组数据突变的功能预测,全外显子组突变的功能预测,剪切变异体的功能预测,疾病相关突变,突变确认等,如千人基因组数据库变异频率进行过滤,各种有害性打分软件打分过滤,各种外显子数据库中的变异频率进行过滤。

4.3.1 使用1000Genomes数据库进行频率注释

下面使用千人基因组数据进行注释:

[kaiwang@biocluster ~/]$ annotate_variation.pl -filter -dbtype 1000g2012apr_eur -buildver hg19 -out ex1 example/ex1.avinput humandb/
NOTICE: Variants matching filtering criteria are written to ex1.hg19_EUR.sites.2012_04_dropped, other variants are written to ex1.hg19_EUR.sites.2012_04_filtered
NOTICE: Processing next batch with 15 unique variants in 15 input lines
NOTICE: Database index loaded. Total number of bins is 2766067 and the number of bins to be scanned is 12
NOTICE: Scanning filter database humandb/hg19_EUR.sites.2012_04.txt...Done

###已存在于数据库中的变异写入 *dropped文件,在数据库中不存在的变异信息将会被写入到*filtered文件.
###第一列为注释数据库名字,第二列为等位基因的突变频率
####需要注意的是,我们也可以使用-maf 0.05 -reverse过滤掉高于0.05的变异;但是过滤ALT等位基因的频率,我们更提倡使用-score_threshold参数。

[kaiwang@biocluster ~/]$ cat ex1.hg19_EUR.sites.2012_04_dropped
1000g2012apr_eur 0.04 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
1000g2012apr_eur 0.87 1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
1000g2012apr_eur 0.81 1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
1000g2012apr_eur 0.06 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
1000g2012apr_eur 0.54 1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
1000g2012apr_eur 0.96 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
1000g2012apr_eur 0.05 16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
1000g2012apr_eur 0.01 16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
1000g2012apr_eur 0.01 16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
1000g2012apr_eur 0.53 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
4.3.2 使用dbsnp数据库进行注释
[kaiwang@biocluster ~/]$ annotate_variation.pl -filter -out ex1 -build hg19 -dbtype snp138 example/ex1.avinput humandb/
NOTICE: Variants matching filtering criteria are written to ex1.hg19_snp138_dropped, other variants are written to ex1.hg19_snp138_filtered
NOTICE: Processing next batch with 15 unique variants in 15 input lines
NOTICE: Database index loaded. Total number of bins is 2858459 and the number of bins to be scanned is 12
NOTICE: Scanning filter database humandb/hg19_snp138.txt...Done

###在dbsnp中已有编号的写入*dropped文件中,在dbsnp中没有的变异写入*filtered文件中

[kaiwang@biocluster ~/]$ cat ex1.hg19_snp138_dropped
snp138 rs35561142 1 11403596 11403596 - AT comments: rs35561142, a 2-bp insertion
snp138 rs149123833 1 1404001 1404001 G T comments: rs149123833, a SNP in 3' UTR of ATAD3C
snp138 rs1000050 1 162736463 162736463 C T comments: rs1000050, a SNP in Illumina SNP arrays
snp138 rs1287637 1 5935162 5935162 A T comments: rs1287637, a splice site variant in NPHP4
snp138 rs11209026 1 67705958 67705958 G A comments: rs11209026 (R381Q), a SNP in IL23R associated with Crohn's disease
snp138 rs6576700 1 84875173 84875173 C T comments: rs6576700 or SNP_A-1780419, a SNP in Affymetrix SNP arrays
snp138 rs15842 1 948921 948921 T C comments: rs15842, a SNP in 5' UTR of ISG15
snp138 rs80338939 13 20763686 20763686 G - comments: rs1801002 (del35G), a frameshift mutation in GJB2, associated with hearing loss
snp138 rs2066844 16 50745926 50745926 C T comments: rs2066844 (R702W), a non-synonymous SNP in NOD2
snp138 rs2066845 16 50756540 50756540 G C comments: rs2066845 (G908R), a non-synonymous SNP in NOD2
snp138 rs2066847 16 50763778 50763778 - C comments: rs2066847 (c.3016_3017insC), a frameshift SNP in NOD2
snp138 rs2241880 2 234183368 234183368 A G comments: rs2241880 (T300A), a SNP in the ATG16L1 associated with Crohn's disease
#*dropped文件
第一列如region-based注释的结果一样以数据库命名;
第二列为已经在数据库的突变的indentifier号;
第三列开始同样是输入文件的内容

参考:
ANNOVAR的使用
ANNOVAR Documentation

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343

推荐阅读更多精彩内容