搜索
查看: 5081|回复: 17

本人的捕获测序数据分析代码,如有错误请大家提出建议

[复制链接]

5

主题

18

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
发表于 2016-11-3 21:33:36 | 显示全部楼层 |阅读模式
本帖最后由 zw579 于 2016-11-26 10:26 编辑

一、原始数据处理
#!/bin/sh
#PBS -N BWA
#PBS -l nodes=8,walltime=10:00:00
#PBS -q ngs
#PBS -o BWA-C18.log
#PBS -e BWA-C18.err   (qsub任务,qstat)(后台运行命令:nohup qsub 路径/文件 &)

#转换原始数据成fastq格式
/home/zz/software/bcl2fastq-1.8.4/bcl2fastq/install/bin/configureBclToFastq.pl --input-dir ~/data160118/160118_SN1033_0316_AC5H67ACXX/Data/Intensities/BaseCalls/ --output-dir ~/data160118/160118_lane6/ --sample-sheet ~/data160118/160118_SN1033_0316_AC5H67ACXX/Data/Intensities/BaseCalls/SampleSheet160118_lane6.csv --fastq-cluster-count 60000000 --no-eamss
cd ~/data160118/160118_lane6
nohup make -j 6


##1 fastq.gz to sam
bwa mem -M -t 8 -R "@RG\tID:C18\tLB:C18\tSM:C18\tPL:ILLUMINA" /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta /home/wcl/data160118/160118_lane6/Project_QSY/Sample_C18/C18_ATCCTG_L006_R1_001.fastq.gz /home/wcl/data160118/160118_lane6/Project_QSY/Sample_C18/C18_ATCCTG_L006_R2_001.fastq.gz |gzip -3 > /home/wcl/ZW/C18/C18.align.sam


首先使用bwa把pair-end双向数据(.fq格式)与参考序列分别进行比对,并把生成的sai文件结合,生成sam文件,然后利用samtools工具把sam文件转换成bam格式,最后利用picard工具按照染色体坐标排列bam文件(Shead表示样本名称)。

-MMark shorter split hits as secondary (for Picard compatibility)
-t INTNumber of threads  线程数

-R STRComplete read group header line. ’\t’ can be used in STR and will be converted to a TAB in the output SAM. The read group ID will be attached to every read in the output. An example is ’@RG\tID:foo\tSM:bar’. [null]

-R int:此参数只应用于pair end中,当没有出现大于此值的最佳比对结果时,将会降低标 准再次进行比对。增加这个值可以提高配对比对的准确率,但是同时会消耗更长的时间,默认是32。




##2 sam to bam
samtools view -@ 8 -bS /home/wcl/ZW/C18/C18.align.sam -o /home/wcl/ZW/C18/C18.align.bam
-@ option can be used to allocate additional threads to be used for compression
-b: indicates that the output is BAM.
-S: indicates that the input is SAM.
-o: specifies the name of the output file.


##3 reorder bam file
java -jar /home/wcl/software/picard-tools-1.138/picard.jar ReorderSam I= /home/wcl/ZW/C18/C18.align.bam O= /home/wcl/ZW/C18/C18.align.reorder.bam R=/home/wcl/annotation/GATK_anno/ucsc.hg19.fasta VALIDATION_STRINGENCY=LENIENT


##4 sort the file
samtools sort /home/wcl/ZW/C18/C18.align.reorder.bam -o /home/wcl/ZW/C18/C18.align.reorder.sorted.bam
sort puts reads without coordinates at the end of the file



##4-1 index the file
samtools index /home/wcl/ZW/C18/C18.align.reorder.sorted.bam


##5 markdup-Sample            注意:本实验外显子捕获用的PCR技术,因此此步骤略去不做,否则reads会被过滤掉很多,导致DP数大大减少
java -jar /home/wcl/software/picard-tools-1.138/picard.jar MarkDuplicates I= /home/wcl/ZW/C18/C18.align.reorder.sorted.bam O= /home/wcl/ZW/C18/C18.align.reorder.sorted.markdup.bam ASSUME_SORTED=true METRICS_FILE= C18.align.metrics VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=true
samtools rmdup does not remove interchromosomal duplicates. MarkDuplicates does remove these duplicates,标记PCR复制



#GATK分析步骤
#Indel Realignment
java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T RealignerTargetCreator \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-I /home/wcl/ZW/C18/C18.align.reorder.sorted.markdup.bam  \
-known /home/wcl/annotation/GATK_anno/1000G_phase1.indels.hg19.vcf \
-known /home/wcl/annotation/GATK_anno/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-o /home/wcl/ZW/C18/C18.realigner.intervals
(通过运行RealignerTargetCreator来确定要进行重新比对的区域)


java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T IndelRealigner \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-I /home/wcl/ZW/C18/C18.align.reorder.sorted.markdup.bam  \
-known /home/wcl/annotation/GATK_anno/1000G_phase1.indels.hg19.vcf \
-known /home/wcl/annotation/GATK_anno/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-targetIntervals /home/wcl/ZW/C18/C18.realigner.intervals \
-o /home/wcl/ZW/C18/C18.realigner.bam

(通过运行IndelRealigner在这些区域内进行重新比对,测序仪产生的质量数据并不总是精确的,存在一些错误,因此需要对这些分、值重新校准,以便获得质量高的变异)


#Base Recalibration(BQSR,碱基质量得分重新校准)

/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T BaseRecalibrator \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-I /home/wcl/ZW/C18/C18.realigner.bam \
-knownSites /home/wcl/annotation/GATK_anno/dbsnp_138.hg19.vcf \
-knownSites /home/wcl/annotation/GATK_anno/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-knownSites /home/wcl/annotation/GATK_anno/1000G_phase1.indels.hg19.vcf \
-o /home/wcl/ZW/C18/output/C18.recal.table

(利用工具BaseRecalibrator,根据一些known sites,生成一个校正质量值所需要的数据文件,GATK网站以“.grp”为后缀命名)




java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T PrintReads \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-I /home/wcl/ZW/C18/C18.realigner.bam  \
-BQSR /home/wcl/ZW/C18/output/C18.recal.table \
-o /home/wcl/ZW/C18/output/C18.recal.bam

(利用工具PrintReads将经过质量值校正的数据输出到新的bam文件中,用于后续的变异检测。)



/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T BaseRecalibrator \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-I /home/wcl/ZW/C18/C18.realigner.bam \
-knownSites /home/wcl/annotation/GATK_anno/dbsnp_138.hg19.vcf \
-knownSites /home/wcl/annotation/GATK_anno/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
-knownSites /home/wcl/annotation/GATK_anno/1000G_phase1.indels.hg19.vcf \
-BQSR /home/wcl/ZW/C18/output/C18.recal.table \
-o /home/wcl/ZW/C18/output/C18.after_recal.table

(这一步主要是为了与校正之前的数据进行比较,最后生成碱基质量值校正前后的比较图,如果不想生成最后BQSR比较图,这一步可以省略)


/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T AnalyzeCovariates \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-before /home/wcl/ZW/C18/output/C18.recal.table \
-after /home/wcl/ZW/C18/output/C18.after_recal.table \
-plots /home/wcl/ZW/C18/output/C18.recalQC.pdf

(这一步会生成评估前后碱基质量值的比较结果,可以选择使用图片和表格的形式展示)


二、 variant    calling
#Germline    variant    calling
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-T HaplotypeCaller \
-I /home/wcl/ZW1/C100/C100.recal.bam \
--dbsnp /home/wcl/annotation/GATK_anno/dbsnp_138.hg19.vcf \
-stand_call_conf 30 \                 #这个是默认值,可以不加
-stand_emit_conf 10 \
-o /home/wcl/ZW1/C100/C100.raw.snps.indels.vcf



#硬质控(hard filters)-将SNP和INDEL分开,然后分别过滤
SNP
#1 Extract the SNPs
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T SelectVariants \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/C100.raw.snps.indels.vcf \
-selectType SNP \
-o /home/wcl/ZW1/C100/output1/C100.raw_snps.vcf


#2 Apply the filter to the SNP call set
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/output1/C100.raw_snps.vcf \
--filterExpression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || HaplotypeScore > 13.0 || QUAL < 30.0 || DP < 5" \
--filterName "my_snp_filter" \
-o /home/wcl/ZW1/C100/output1/C100.filtered_snps.vcf

#去掉没通过质控的位点
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T SelectVariants \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/output1/C100.filtered_snps.vcf \
--excludeFiltered \
-o /home/wcl/ZW1/C100/output1/C100.filtered_PASS.snps.vcf
#评估SNP
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \               
-T VariantEval \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
--dbsnp /home/wcl/annotation/GATK_anno/dbsnp_138.hg19.vcf \
-eval /home/wcl/ZW1/C100/output1/C100.filtered_PASS.snps.vcf \
-o /home/wcl/ZW1/C100/output1/C100.filtered_PASS.snps.vcf.eval




INDEL
#Extract the INDELS
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T SelectVariants \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/C100.raw.snps.indels.vcf \
-selectType INDEL \
-o /home/wcl/ZW1/C100/output2/C100.raw_indels.vcf


#Apply the filter to the Indel call set
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T VariantFiltration \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/output2/C100.raw_indels.vcf \
--filterExpression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0 || QUAL < 30.0 || DP < 5" \
--filterName "my_indel_filter" \
-o /home/wcl/ZW1/C100/output2/C100.filtered_indels.vcf


/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T SelectVariants \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW1/C100/output2/C100.filtered_indels.vcf \
--excludeFiltered \
-o /home/wcl/ZW1/C100/output2/C100.filtered_PASS.indels.vcf




三、变异检测
由于上一步将SNP与INDEL分开,所以这一步分开检测,生成两份csv文件 (table_annovar.pl可以输入VCF文件不用做格式转换这步,但是会输出同名的vcf文件,可以加参数改变)
{ SNP
#Format Conversion
/home/wcl/software/annovar/convert2annovar.pl -format vcf4 /home/wcl/ZW1/C100/output1/C100.filtered_PASS.snps.vcf -outfile /home/wcl/ZW1/C100/output1/C100.snps.avinput -includeinfo



#annotation
/home/wcl/software/annovar/table_annovar.pl \
/home/wcl/ZW1/C100/output1/C100.snps.avinput \
/home/wcl/software/annovar/humandb/ -buildver hg19 \
-out /home/wcl/ZW1/C100/output2/C100snp -remove \
-protocol refGene,cytoBand,exac03,esp6500siv2_all,clinvar_20140929,1000g2014oct_all,snp138,ljb26_all \
-operation g,r,f,f,f,f,f,f \
-nastring NA -csvout -otherinfo



INDEL
#Format Conversion

/home/wcl/software/annovar/convert2annovar.pl -format vcf4 /home/wcl/ZW1/C100/output2/C100.filtered_PASS.indels.vcf \
-outfile /home/wcl/ZW1/C100/output2/C100.indels.avinput -includeinfo



#annotation


/home/wcl/software/annovar/table_annovar.pl \
/home/wcl/ZW1/C100/output2/C100.indels.avinput \
/home/wcl/software/annovar/humandb/ -buildver hg19 \
-out /home/wcl/ZW1/C100/output2/C100indel -remove \
-protocol refGene,cytoBand,exac03,esp6500siv2_all,clinvar_20140929,1000g2014oct_all,snp138,ljb26_all \
-operation g,r,f,f,f,f,f,f \
-nastring NA -csvout -otherinfo

}
#推荐用VCF文件输入和输出注释,方便后面导入plink分析
/home/wcl/software/annovar/table_annovar.pl \
/home/wcl/ZW1/Q50/output/Q50.filtered_PASS.snps.vcf \
/home/wcl/software/annovar/humandb/ -buildver hg19 \
-out /home/wcl/ZW/csv/Q50snp -remove \
-protocol refGene,cytoBand,1000g2014oct_all,snp138 \
-operation g,r,f,f \
-nastring . -vcfinput -otherinfo


#合并VCF进行后续plink分析
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T CombineVariants \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \

--variant /home/wcl/ZW/csv/Q8snp.hg19_multianno.vcf \
--variant /home/wcl/ZW/csv/Q9snp.hg19_multianno.vcf \
-o /home/wcl/ZW/csv/output/put.vcf \
-genotypeMergeOptions UNIQUIFY



#转换vcf成bim文件之前对vcf文件进行整理
sed -i "/chrM/d" /home/wcl/ZW/csv/output/put.vcf     #删除chrM所在的行
sed -i "s/chr//g" /home/wcl/ZW/csv/output/output.bim #将chr换成空




#chrN_random和chrUn:基因组文件中通常含有类似chr9_gl000198_random和chrUn_gl000211的基因组。根据UCSC解释,chrN_random包括基因组已知但具体位置未知序列,或者位置已知但具体内容未完成序列。ChrUn中包括一些具体位置未知的序列。
sed -i "/chrM/d" /home/wcl/ZW/csv/output/put.vcf \

#用VariantsToBinaryPed工具将VCF转换成plink所用格式                                                                                                                     
#注意:输入的metadata样本名称和顺序必须和vcf文件一样,否则不能识别,切记切记!
/home/wcl/software/jdk1.7.0_80/bin/java -jar /home/wcl/software/GATK/GATK-3.4.46/GenomeAnalysisTK.jar \
-T VariantsToBinaryPed \
-R /home/wcl/annotation/GATK_anno/ucsc.hg19.fasta \
-V /home/wcl/ZW/csv/output/put.vcf \
-m /home/wcl/ZW/csv/output/metadata.fam \
-bed /home/wcl/ZW/csv/output/output.bed \
-bim /home/wcl/ZW/csv/output/output.bim \
-fam /home/wcl/ZW/csv/output/output.fam \
--minGenotypeQuality 20


#转格式之后要把bim文件中的chr换掉,否则plink不能识别文件

PLINK使用:
#将二进制bed文件转换成ped文件
./plink --bfilegwas_example --recode --out gwas_example

#对样本和SNP个数进行计数
wc -l gwas_example.bim            #本次奥氮平所有样本总的snp个数:520507 SNPs
#将bim文件中的chr换成空,否则会报错
sed -i "s/chr//g" /home/wcl/ZW/csv/output/output.bim
wc -l gwas_example.fam           #样本数:78

#样本质控

/home/wcl/software/plink-1.07-x86_64/plink --bfile /home/wcl/ZW/csv/output/output --missing --out /home/wcl/ZW/csv/output/sample

#此步骤会生成gwas_example.lmiss gwas_example.imiss两个文件。   gwas_example.imiss lists a sample per row with genotyping call rate information
less -S gwas_example.imiss
file:///E:/%E6%9C%89%E9%81%93%E4%BA%91%E7%AC%94%E8%AE%B0/zhouweihkd@163.com/d36f6a1d86d2499098b1acaaa6627bd2/clipboard.png

#How many samples have genotype missing rate > 0.02?   
awk '$6 > 0.02' gwas_example.imiss


#SNP位点质控
#去掉missing genotype rate >0.05的位点
awk '$5 > 0.05' /home/wcl/ZW/csv/output/sample.lmiss
awk 'NR < 2 { next } $5 > 0.05' /home/wcl/ZW/csv/output/sample.lmiss > /home/wcl/ZW/csv/output/sample.lmiss.exclude
/home/wcl/software/plink-1.07-x86_64/plink --bfile /home/wcl/ZW/csv/output/sample --exclude /home/wcl/ZW/csv/output/sample.lmiss.exclude --make-bed --out /home/wcl/ZW/csv/output/sample-qc2


[wcl@console output]$wc -l output2qc1.bim
908 output2qc1.bim
[wcl@console output]$wc -l output2.lmiss
17981 output2.lmiss
[wcl@console output]$wc -l output2.lmiss.exclude
17072 output2.lmiss.exclude
#HWE QC
/home/wcl/software/plink-1.07-x86_64/plink --bfile /home/wcl/ZW/csv/output/output2qc1 --hardy --out /home/wcl/ZW/csv/output/output2qc1
awk '$3=="UNAFF" && $9 < 0.001 {print $2}' /home/wcl/ZW/csv/output/output2qc1.hwe > /home/wcl/ZW/csv/output/output2qc1.hwe.exclude
/home/wcl/software/plink-1.07-x86_64/plink --bfile /home/wcl/ZW/csv/output/output2qc1 --exclude /home/wcl/ZW/csv/output/output2qc1.hwe.exclude --make-bed --out /home/wcl/ZW/csv/output/output2qc2


#协变量分析
‘plink --file clean -baihui14 --logistic --covar name.covar --out name
cova.xlsx 改成name.covar


#关联分析
/home/wcl/software/plink-1.07-x86_64/plink --file /home/wcl/ZW/csv/output/output2qc2 --assoc --counts --ci 0.95 --out /home/wcl/ZW/csv/output/output2assoc






上一篇:增强子在特定的细胞系状态下是否是激活的?
下一篇:【学分析】宏基因组分析教程分享
回复

使用道具 举报

5

主题

18

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
 楼主| 发表于 2016-11-8 13:50:23 | 显示全部楼层
修正:plink命令加入--allow-extra-chr 就可以识别了。
回复 支持 反对

使用道具 举报

11

主题

50

帖子

275

积分

版主

Rank: 7Rank: 7Rank: 7

积分
275
发表于 2016-11-11 04:39:58 | 显示全部楼层
好长,有些指令写个循环吧,赞分享!
回复 支持 反对

使用道具 举报

1

主题

7

帖子

90

积分

注册会员

Rank: 2

积分
90
发表于 2016-11-17 08:58:05 | 显示全部楼层
外显子捕获用PCR测序 得到的数据不用picard去重么?
回复 支持 反对

使用道具 举报

5

主题

18

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
 楼主| 发表于 2016-11-17 13:44:04 | 显示全部楼层
Mint 发表于 2016-11-11 04:39
好长,有些指令写个循环吧,赞分享!

嗯,还在学习中,有些不会
回复 支持 反对

使用道具 举报

5

主题

18

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
 楼主| 发表于 2016-11-17 13:45:55 | 显示全部楼层
kkshaxqd 发表于 2016-11-17 08:58
外显子捕获用PCR测序 得到的数据不用picard去重么?

去重后发现测序深度大大下降,所以就没去了,这也是我不明白的地方,请指教。
回复 支持 反对

使用道具 举报

1

主题

7

帖子

90

积分

注册会员

Rank: 2

积分
90
发表于 2016-11-19 09:31:21 | 显示全部楼层
zw579 发表于 2016-11-17 13:45
去重后发现测序深度大大下降,所以就没去了,这也是我不明白的地方,请指教。 ...

应该还是要去的,这个属于干扰吧。不过有时去的的确很多,从1万多reads直接到几千。
回复 支持 反对

使用道具 举报

0

主题

3

帖子

60

积分

注册会员

Rank: 2

积分
60
发表于 2016-11-21 13:44:22 | 显示全部楼层
请教一个问题,你这是用的annovar对vcf进行注释的对吧? annovar 需要先用convert2annovar.pl 把vcf文件转化成avinput,注释完毕后是需要用什么程序再转化回vcf文件呢?不太理解,谢谢!
回复 支持 反对

使用道具 举报

5

主题

18

帖子

272

积分

中级会员

Rank: 3Rank: 3

积分
272
 楼主| 发表于 2016-11-26 10:23:29 | 显示全部楼层
皖城小西 发表于 2016-11-21 13:44
请教一个问题,你这是用的annovar对vcf进行注释的对吧? annovar 需要先用convert2annovar.pl 把vcf文件转 ...

#推荐用VCF文件输入和输出注释,方便后面导入plink分析
/home/wcl/software/annovar/table_annovar.pl \
/home/wcl/ZW1/Q50/output/Q50.filtered_PASS.snps.vcf \
/home/wcl/software/annovar/humandb/ -buildver hg19 \
-out /home/wcl/ZW/csv/Q50snp -remove \
-protocol refGene,cytoBand,1000g2014oct_all,snp138 \
-operation g,r,f,f \
-nastring . -vcfinput -otherinfo
这一步里输出的是vcf文件吧
回复 支持 反对

使用道具 举报

1

主题

3

帖子

29

积分

新手上路

Rank: 1

积分
29
发表于 2016-12-29 09:56:11 | 显示全部楼层
赞。不过我的去重了,丢了20%的reads。。。
回复 支持 反对

使用道具 举报

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

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-8-16 01:01 , Processed in 0.124045 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.