在1000genome的FTP服务器上可以下载一个all的vcf文件，里面可以看到AFR, AMR, EAS, EUR, SAS人群的allele频率，但是该种族下面的亚群的频率信息需要在http://grch37.ensembl.org搜索得到，比如 http://grch37.ensembl.org/Homo_sapiens/Variation/Population?db=core;r=1:230845294-230846294;v=rs699;vdb=variation;vf=102788013 ，还有一种方式，就是下载包含所有样本的突变信息的VCF文件，利用vcftools计算。

The allele frequency of super population (AFR, AMR, EAS, EUR, SAS, see http://www.1000genomes.org/category/population/) can be obtained from all.vcf.

However, the allele population frequency in subpopulation is not well obtained.

One way is search the web http://grch37.ensembl.org by rs identifier, E.g. http://grch37.ensembl.org/Homo_sapiens/Variation/Population?db=core;r=1:230845294-230846294;v=rs699;vdb=variation;vf=102788013 .

Another way is calculated on your local machine. The following will introduce how to get allele frequency of CHB population in chr1 chromosome. CHB Han Chinese in Beijing

修改MySQL远程访问权限

1

grant all privileges on *.* to 'root'@'%' identified by '123456';

# . 表示允许所有表 ‘root’ 表示允许root用户 ‘%’ 表示任意ip ‘123456’ 表示密码

1

flush privileges;

刷新

bismark调用bowtie2进行比对，调用samtools生成bam文件，因此在运行bismark之前，需要安装bowtie2和samtools

请注意，fastq文件要进行质控，比如去掉低质量的reads，去掉adaptor等，可以看本文最下方推荐的PPT,本文不介绍，此外本本只介绍到序列比对，后续的统计分析没有介绍，有兴趣的朋友可以关注swDMR和methykit工具包。

安装bismark

1
2

wget http://www.bioinformatics.babraham.ac.uk/projects/bismark/bismark_v0.16.1.tar.gz
tar -xvzf bismark_v0.16.1.tar.gz

生成转换后的基因组

# --path\_to\_bowtie后面跟的是文件夹
# --verbose 输出log信息
# ./ref 文件夹中有一个基因组fasta文件
# --bowtie2指明用的是bowtie2
./bismark\_v0.16.1/bismark\_genome\_preparation --path\_to_bowtie /home/zzx/bowtie2-2.2.9/ --bowtie2 --verbose ./ref/

因为在重亚硫酸盐甲基化测序中，因为未甲基化的C会变为T，在正链表现为C–>T，但是在负链有C变为T，转换为正链时，即为G–>A，所以基因组需要进行两种转化，才能用于比对。在基因组目录下产生Bisulfite_Genome目录，有CT_conversion和GA_conversion文件夹，这两个文件夹包含转换后的fasta文件和bowtie2建立的索引bt2文件。

fastq中的BS转换后的read与转换的参考基因组比对，得到在参考基因组中的位置，再与原始的参考基因组比较，确定methylate call

运行trimmomatic的默认参数

1

java -jar trimmomatic.jar PE -threads 16 -phred33 sample_1.fastq sample_2.fastq sample_trimmed_paired_1.fastq.gz sample_trimmed_unpaired_1.fastq.gz sample_trimmed_paired_2.fastq.gz sample_trimmed_unpaired_2.fastq.gz ILLUMINACLIP:adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36

输出文件中

sample_trimmed_unpaired_1.fastq.gz是sample_trimmed_unpaired_2.fastq.gz的十多倍，异常的大。1文件（forward reads）中unpaired reads非常多，显著多余2文件（reverse reads）中的unpaired reads。