参考基因组的准备
-
UCSC Genome Browser
主页 > Downloads > Genome Data > Human
上面两个都可以
点进去在页面底部就可以看到hg38参考基因组了。
#比如这个
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz
-
Ensembl
界面简洁美观
点进去之后选择Download DNA sequence就可以看到下载界面了。
#依次下载22条常染色体和1对性染色体,然后cat合并就可以了
wget ftp://ftp.ensembl.org/pub/release-94/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.chromosome.1.fa.gz
-
GATK官网(我用的这个)
主页 > Download > Resource Bundle在页面最下方可以看到FTP Server Access, 点进去可以看到:
链接:ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/,在上级目录中也可以找到hg19的参考文件
#dict fai两个索引文件在后面局部重比对目标创建的过程中会用到,没有会这样报错
##### ERROR MESSAGE: Fasta index file XXX.fa.fai for reference XXX.fa does not exist.
##### Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.
##### ERROR MESSAGE: Fasta dict file XXX.dict for reference XXX.fa does not exist.
##### Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.
#另外这两个命令是可以生成dict fai两个索引文件的,忘了的话可以补上
samtools faidx XXX.fa
java -jar picard.jar CreateSequenceDictionary R=XXX.fa
#还有几个vcf文件也需要下载,后面的局部重比对,碱基校正,变异检测,变异校正都会以此为参考。
1000G_phase1.snps.high_confidence.hg38.vcf.gz
1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
dbsnp_146.hg38.vcf.gz
dbsnp_146.hg38.vcf.gz.tbi
hapmap_3.3.hg38.vcf.gz
hapmap_3.3.hg38.vcf.gz.tbi
Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi
为参考基因组建立索引
一本书有了目录之后可以更快地查找书中的内容,参考基因组的索引也是这个作用。
#简单查看一下参考基因组的信息
$ grep ">" Homo_sapiens_assembly38.fasta | less -S -x 16
>chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38
>chr2 AC:CM000664.2 gi:568336022 LN:242193529 rl:Chromosome M5:f98db672eb0993dcfdabafe2a882905c AS:GRCh38
...
#创建索引
bwa index -a bwtsw Homo_sapiens_assembly38.fasta
#会生成amb、ann、bwt、pac、sa结尾的5个索引文件文件
#参考基因组和索引文件放在一个文件夹里
构建索引时需要注意的问题:bwa构建索引有两种算法,两种算法都是基于BWT的,这两种算法通过参数-a is 和-a bwtsw进行选择。其中-a bwtsw对于短的参考序列是不工作的,必须要大于等于10Mb;-a is是默认参数,这个参数不适用于大的参考序列,必须要小于等于2G。
网友评论