200826 Circ之旅3-构建人类基因组索引
注：后面可能还会构建鼠源的x
参考：
bwa bowtie2 salmon subread hisat2建索引和比对
「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
步骤和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
但是问题在于解压出来的文件，大小是在有点奇怪啊。
-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
最后是在服务器上跑完了。所以服务器还是很重要的，自家电脑大部分时候一是可能会跑一半宕机，还有可能是撑满内存。
