200826 Circ之旅3-构建人类基因组索引

注：后面可能还会构建鼠源的x
参考：

RNA-seq(5):序列比对：Hisat2

人类基因组hg19、hg38构建bwa索引

生物信息学习——bowtie2使用手册

bwa bowtie2 salmon subread hisat2建索引和比对

史上最快的转录组流程-subread

生信软件 | bowtie2（测序序列与参考序列比对）

「RNA-seq分析软件」RNA-seq比对工具STAR学习笔记

我这用的hg19和hg38，构建了几个软件的，bwa, hisat2, subread,bowtie2,STAR
STAR建议在服务器上跑

获取测序源文件

hg19

hg19的UCSC下载链接
最后反正我直接下了总的fa文件，也有让下别的。我偷懒
下以下几个文件：

hg19.fa.gz - "Soft-masked" assembly sequence in one file. Repeats from RepeatMasker and Tandem Repeats Finder (with period of 12 or less) are shown in lower case; non-repeating sequence is shown in upper case.
md5sum.txt - checksums of files in this directory

检测文件完整性：感觉有点类似哈希校验那种感觉。

cat md5sum.txt
# 将hg19.fa.gz和对应的校验码提出来
echo "806c02398f5ac5da8ffd6da2d1d5d1a9  hg19.fa.gz" > check_md5sum.txt
# 用md5sum校验
md5sum -c check_md5sum.txt
# 如果校验完成，显示 hg19.fa.gz: OK

注意这里面写的是>，不是>>。好像>是覆盖写入，>>是在文件末尾加。

校验没问题就解压

gzip -dk hg19.fa.gz
-d, --decompress  decompress #不加d就是压缩了
-k, --keep        keep (don not delete) input files #保存源文件存在，不加源文件会被删除。

解压之后获得hg19.fa然后对这个文件建立bwa的索引就行。

hg38

hg38的UCSC下载链接

步骤和hg19完全一样，然后下载以下文件

hg38.fa.gz - "Soft-masked" assembly sequence in one file. Repeats from RepeatMasker and Tandem Repeats Finder (with period of 12 or less) are shown in lower case; non-repeating sequence is shown in upper case.
md5sum.txt - checksums of files in this directory

gtf文件

hg38的各种gtf文件
hg19的各种gtf文件
可能可以下ens或者know，反正我下的knowGene

bwa

hg19

语法：

bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>
bwa index -a bwtsw hg19.fa
-p STR   输出数据库的前缀；默认和输入的文件名一致，输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀。
-a [is|bwtsw]   构建index的算法，有两个算法： is 是默认的算法，虽然相对较快，但是需要较大的内存，当构建的数据库大于2GB的时候就不能正常工作了。 bwtsw 对于短的参考序列式不工作的，必须要大于等于10MB, 但能用于较大的基因组数据，比如人的全基因组。

输结果如下：

[BWTIncConstructFromPacked] 680 iterations done. 6241341392 characters processed.
[BWTIncConstructFromPacked] 690 iterations done. 6264217232 characters processed.
[bwt_gen] Finished constructing BWT in 695 iterations.
[bwa_index] 1878.26 seconds elapse.
[bwa_index] Update BWT... 15.69 sec
[bwa_index] Pack forward-only FASTA... 15.37 sec
[bwa_index] Construct SA from BWT and Occ... 941.87 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index -a bwtsw hg19.fa
[main] Real time: 2905.553 sec; CPU: 2872.593 sec

最后会生成五个文件：hg19.fa.amb、hg19.fa.ann、hg19.fa.bwt、hg19.fa.pac、hg19.fa.sa

具体这五个文件是啥啥啥就不管了，留个参考链接：BWA源码阅读笔记（二）索引文件amb/ann/pac文件是什么?

hg38

步骤和hg19完全一样

bwa index -a bwtsw hg38.fa

hisat2

这个我是直接去官网下的，链接：hisat2索引下载地址

对应的是UCSC hg19和UCSC hg38，至于和剩下两个GRCh38和GRCh37的区别，见下面连接：

hg19、GRCH37、b37、hs37d5介绍和区别

hg19

但是问题在于解压出来的文件，大小是在有点奇怪啊。

-rw-r--r-- 1 dick dick 926M 11月 20  2015 genome.1.ht2
-rw-r--r-- 1 dick dick 691M 11月 20  2015 genome.2.ht2
-rw-r--r-- 1 dick dick 4.8K 11月 20  2015 genome.3.ht2
-rw-r--r-- 1 dick dick 691M 11月 20  2015 genome.4.ht2
-rw-r--r-- 1 dick dick 1.2G 11月 20  2015 genome.5.ht2
-rw-r--r-- 1 dick dick 704M 11月 20  2015 genome.6.ht2
-rw-r--r-- 1 dick dick    8 11月 20  2015 genome.7.ht2
-rw-r--r-- 1 dick dick    8 11月 20  2015 genome.8.ht2
-rwxr-xr-x 1 dick dick 1.3K 11月 20  2015 make_hg19.sh

这个大小就不是非常让人信服，所以决定用我自己的hisat重新跑一下。

hisat2-build [options]* <reference_in> <ht2_index_base>
reference_in            comma-separated list of files with ref sequences
hisat2_index_base       write ht2 data to files with this dir/basename
-p <int>                number of threads
hisat2-build -p 16 hg19.fa genome

最后输出：

-rw-r--r-- 1 dick dick 926M 6月  22 17:05 genome.1.ht2
-rw-r--r-- 1 dick dick 691M 6月  22 17:05 genome.2.ht2
-rw-r--r-- 1 dick dick 4.8K 6月  22 16:52 genome.3.ht2
-rw-r--r-- 1 dick dick 691M 6月  22 16:52 genome.4.ht2
-rw-r--r-- 1 dick dick 1.2G 6月  22 17:07 genome.5.ht2
-rw-r--r-- 1 dick dick 704M 6月  22 17:07 genome.6.ht2
-rw-r--r-- 1 dick dick   12 6月  22 16:52 genome.7.ht2
-rw-r--r-- 1 dick dick    8 6月  22 16:52 genome.8.ht2

完全一样，因该就是这样一个索引了。

hg38

基本上一样，就记录一下最后的输出吧。

-rw-r--r-- 1 dick 974M 11月 20  2015 genome.1.ht2
-rw-r--r-- 1 dick 728M 11月 20  2015 genome.2.ht2
-rw-r--r-- 1 dick  15K 11月 20  2015 genome.3.ht2
-rw-r--r-- 1 dick 728M 11月 20  2015 genome.4.ht2
-rw-r--r-- 1 dick 1.3G 11月 20  2015 genome.5.ht2
-rw-r--r-- 1 dick 741M 11月 20  2015 genome.6.ht2
-rw-r--r-- 1 dick    8 11月 20  2015 genome.7.ht2
-rw-r--r-- 1 dick    8 11月 20  2015 genome.8.ht2
-rwxr-xr-x 1 dick 1.3K 11月 20  2015 make_hg38.sh

Total time for call to driver() for forward index: 00:16:53

subread

这个软件好像还比较简单，一句就可以了。

subread-buildindex -o hg19 hg19.fa

这个程序写的不是很好，-o后面更的是要生成index的名字，没搞懂 一开始以为是输出文件夹，然后文档写的也不清楚。

生成文件如下：

-rw-r--r-- 1 dick dick 749M 6月  24 22:59 hg19.00.b.array
-rw-r--r-- 1 dick dick 4.9G 6月  24 22:59 hg19.00.b.tab
-rw-r--r-- 1 dick dick 5.5K 6月  24 22:57 hg19.files
-rw-r--r-- 1 dick dick    0 6月  24 22:46 hg19.log
-rw-r--r-- 1 dick dick 2.3K 6月  24 22:59 hg19.reads

hg38也是一样的。

ubread-buildindex -o hg38 hg38.fa

bowtie2

网上似乎可以下载：Bowtie 2: Manual 下载在右边的栏目里，文档好像这儿也可以找到，这个软件的文档写的非常可以。

但是只有hg19和GRCh38，所以我用的hg38，就手动建立了。

命令行如下：

bowtie2-build --threads 16 hg19.fa hg19
bowtie2-build --threads 16 hg38.fa hg38
这条命令在结束前应该会打印很多行输出。当其运行完毕时，当前文件夹会产生6个新的文件，它们的文件名都以hg19开始，分别以.1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2和.rev.2.bt2结束。这些文件构成了索引——你完成了！
--thread 10 这个参数是指定多线程的，肯定要用的，其他没研究。（我这破电脑，最多因该可以跑到16）
注意：这个命令是有参数的，在文档里可以找到，懒得去研究了。

绝了，我做出来的索引居然还大一点点。不管了还是先用我自己的索引吧。

star

注意：STAR非常吃内存，所以可以用一下spares mode，我最后是在服务器上完成的。

代码：

STAR  --runMode genomeGenerate \
    --genomeDir ref \
    --runThreadN 20 \
    --genomeFastaFiles reference.fa\
    --sjdbGTFfile reference.gtf

常用参数说明

• --runThreadN 线程数 :设置线程数
• --runMode genomeGenerate : 设置模式为构建索引
• --genomedDir 索引文件存放路径 : 必须先创建文件夹
• --genomeFastaFiles 基因组fasta文件路径 : 支持多文件路径
• --sjdbGTFfile gtf文件路径 : 可选项，高度推荐,用于提高比对精确性
• --sjdbOverhang 读段长度: 后续回帖读段的长度, 如果读长是PE 100， 则该值设为100-1=99

坑1：如果设置--sjdbGTFfile或--sjdbFileChrStartEnd，就需要设置--sjdbOverhang，这个选项参数是测序长度-1好像，印象里我的测序是50的读长，设置49看看。

坑2：这个软件会非常吃内存，如果不限制内存使用，索引可能会使用30G往上的内存然后在sorting Suffix Array chunks and saving them to disk这一步会报错killed，没有任何提示，所以要限制一下内存。

试一下限制ram --limitGenomeGenerateRAM 10000000000，这个命令好蠢啊，是限制多少byte，还要换算，我电脑是16g，那么换算下来大概17000000000byte，这样可以了，限制之后好像速度会变得巨巨巨巨慢。

坑3：限制内存一定是限制到空闲内存的量以下，否则会出现terminate called after throwing an inatance of 'std::bas_alloc'查了一下解决方法好像是用--genomeChrBinNbits = 11

The number here depends on the number of sacffolds/sequences and assembly size and is calculated by using the formula given in the manual asmin(18,log2[max(GenomeLength/NumberOfReferences,ReadLength)])

注：贴个教程：基因组注释文件(GFF,GTF)下载的四种方法

最后是去UCSC下载的，教程里那个下下来不是很好用，建议直接去各个基因的下载页面下载东西。

hg38下载页面里面有个genes文件夹、hg19下载页面里面有个genes文件夹

有ens、known、ncbi、ref四种gtf文件，我下的knowngene，我理解是已知基因序列的GTF。

hg19

STAR  --runMode genomeGenerate \
    --limitGenomeGenerateRAM 10000000000 \
    --genomeDir ~/Circ/index/STAR/hg19/ \
    --genomeFastaFiles hg19.fa\
    --sjdbGTFfile ./gtf/hg19.knownGene.gtf\
    --sjdbOverhang 49 --runThreadN 16

记录一下输出：

流程：
这个软件是先生成很多分段的SA文件，然后所有SA文件再合成一个整的索引
Jul 06 09:39:06 ..... started STAR run
Jul 06 09:39:06 ... starting to generate Genome files
Jul 06 09:40:22 ..... processing annotations GTF
Jul 06 09:40:47 ... starting to sort Suffix Array. This may take a long time...
Jul 06 09:41:07 ... sorting Suffix Array chunks and saving them to disk...
        在这一步非常慢，会生成很多分段的SA文件。
Jul 06 09:58:03 ... loading chunks from disk, packing SA...
Jul 06 10:00:27 ... finished generating suffix array
Jul 06 10:00:27 ... generating Suffix Array index
Jul 06 10:05:45 ... completed Suffix Array index
Jul 06 10:05:45 ..... inserting junctions into the genome indices
Jul 06 10:08:27 ... writing Genome to disk ...
Jul 06 10:08:30 ... writing Suffix Array to disk ...
Jul 06 10:10:51 ... writing SAindex to disk
Jul 06 10:11:02 ..... finished successfully

输出文件

-rw-rw-r-- 1 dick dick  688 7月   6 09:40 chrLength.txt
-rw-rw-r-- 1 dick dick 2.0K 7月   6 09:40 chrNameLength.txt
-rw-rw-r-- 1 dick dick 1.3K 7月   6 09:40 chrName.txt
-rw-rw-r-- 1 dick dick 1021 7月   6 09:40 chrStart.txt
-rw-rw-r-- 1 dick dick  25M 7月   6 09:40 exonGeTrInfo.tab
-rw-rw-r-- 1 dick dick  11M 7月   6 09:40 exonInfo.tab
-rw-rw-r-- 1 dick dick 2.2M 7月   6 09:40 geneInfo.tab
-rw-rw-r-- 1 dick dick 3.0G 7月   6 10:08 Genome
-rw-rw-r-- 1 dick dick  640 7月   6 10:08 genomeParameters.txt
-rw-rw-r-- 1 dick dick  21K 7月   6 10:11 Log.out
-rw-rw-r-- 1 dick dick  23G 7月   6 10:10 SA
-rw-rw-r-- 1 dick dick 1.5G 7月   6 10:10 SAindex
-rw-rw-r-- 1 dick dick 7.4M 7月   6 10:05 sjdbInfo.txt
-rw-rw-r-- 1 dick dick 9.7M 7月   6 09:40 sjdbList.fromGTF.out.tab
-rw-rw-r-- 1 dick dick 6.6M 7月   6 10:05 sjdbList.out.tab
-rw-rw-r-- 1 dick dick 4.8M 7月   6 09:40 transcriptInfo.tab

hg38

STAR  --runMode genomeGenerate \
    --genomeDir ~/Circ/index/STAR/hg38/ \
    --genomeFastaFiles hg38.fa\
    --sjdbGTFfile ./gtf/hg38.knownGene.gtf \
    --sjdbOverhang 49 --runThreadN 25

输出文件：

-rw-rw-r-- 1 dick dick 3.0K 7月   6 10:24 chrLength.txt
-rw-rw-r-- 1 dick dick  12K 7月   6 10:24 chrNameLength.txt
-rw-rw-r-- 1 dick dick 8.5K 7月   6 10:24 chrName.txt
-rw-rw-r-- 1 dick dick 4.9K 7月   6 10:24 chrStart.txt
-rw-rw-r-- 1 dick dick  51M 7月   6 10:24 exonGeTrInfo.tab
-rw-rw-r-- 1 dick dick  21M 7月   6 10:24 exonInfo.tab
-rw-rw-r-- 1 dick dick 9.1M 7月   6 10:24 geneInfo.tab
-rw-rw-r-- 1 dick dick 3.1G 7月   6 11:08 Genome
-rw-rw-r-- 1 dick dick  640 7月   6 11:08 genomeParameters.txt
-rw-rw-r-- 1 dick dick 6.6M 7月   6 11:10 Log.out
-rw-rw-r-- 1 dick dick  24G 7月   6 11:09 SA
-rw-rw-r-- 1 dick dick 1.5G 7月   6 11:10 SAindex
-rw-rw-r-- 1 dick dick  12M 7月   6 11:05 sjdbInfo.txt
-rw-rw-r-- 1 dick dick  17M 7月   6 10:24 sjdbList.fromGTF.out.tab
-rw-rw-r-- 1 dick dick  11M 7月   6 11:05 sjdbList.out.tab
-rw-rw-r-- 1 dick dick  16M 7月   6 10:24 transcriptInfo.tab

最后是在服务器上跑完了。所以服务器还是很重要的，自家电脑大部分时候一是可能会跑一半宕机，还有可能是撑满内存。