HLA-VBseq 利用全基因组测序的数据,可以提供8位的HLA分型结果,其文献链接如下
https://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-16-S2-S7
在该文献中,利用30X的全基因组数据,对HLA-VBSeq, PHLAT, HLAminer这3款软件的分型结果进行了评估,准确率汇总如下
可以看到,只有HLA-VBSeq提供了8位的分型结果,准确率高达99.94%;对于2位到4位的分型结果,其准确率也高于另外两款软件。
同时还评估了不同测序量时,各种软件提供的4位分型结果的准确率,结果如下
在不同条件下,HLA-VBseq的准确率都是最高的。由此可见,该软件的分型效果还是相当不错的,官网如下
http://nagasakilab.csml.org/hla/
该软件采用java语言开发,直接下载HLAVBseq.jar
就可以了,除了该文件之外,还需要下载以下几个文件
bamNameIndex.jar
SamToFastq.jar
parse_result.pl
hla_all.fasta
Allelelist.txt
前三个程序在处理fastq文件时会用到;后两个文件是从IMGA/HLA数据库下载的,如果觉得官网提供的版本较老,可以从IMGA/HLA数据库下载最新版。
软件的步骤较多,首先将fastq序列与参考基因组进行比对,得到bam文件,然后对该bam文件进行操作。步骤如下:
利用samtools view 命令挑选出比对到HLA区域的reads , 命令如下
samtools view -hb align.bam chr6:29907037-29915661 chr6:31319649-31326989 chr6:31234526-31241863 chr6:32914391-32922899 chr6:32900406-32910847 chr6:32969960-32979389 chr6:32778540-32786825 chr6:33030346-33050555 chr6:33041703-33059473 chr6:32603183-32613429 chr6:32707163-32716664 chr6:32625241-32636466 chr6:32721875-32733330 chr6:32405619-32414826 chr6:32544547-32559613 chr6:32518778-32554154 chr6:32483154-32559613 chr6:30455183-30463982 chr6:29689117-29699106 chr6:29792756-29800899 chr6:29793613-29978954 chr6:29855105-29979733 chr6:29892236-29899009 chr6:30225339-30236728 chr6:31369356-31385092 chr6:31460658-31480901 chr6:29766192-29772202 chr6:32810986-32823755 chr6:32779544-32808599 chr6:29756731-29767588 | samtools fastq - -1 R1.fq -2 R2.fq
需要注意的是,在使用view命令时,虽然也可以直接提供一个bed格式的文件来挑选特定区域的reads,但是这种用法不会利用到bam文件的索引,所以速度很慢。对于全基因组数据,bam文件很大,上述写法虽然冗长,但是执行效率高。
利用samtools view 命令挑选出没有比对上参考基因组的reads, 命令如下:
samtools view -hb -f 12 /home/pub/output/WGS/18B0315D/6343/6343_final.bam | samtools fastq - -1 unmapped_R1.fq -2 unmapped_R2.fq
将比对到HLA区域的reads和没比对上参考基因组的reads合并,命令如下
cat R1.fq unmapped_R1.fq > R1.fastqcat R2.fq unmapped_R2.fq > R2.fastq
利用bwa软件,将上一步得到的reads与HLA参考序列比对,命令如下
bwa index hla_all.fastabwa mem -t 8 -P -L 10000 -a hla_all.fasta R1.fastq R2.fastq > out.sam
HLA-VBSeq支持双端或者单端测序的数据,这里以双端数据为例,用法如下
java -jar HLAVBSeq.jar hla_all.fasta out.sam result.txt --alpha_zero 0.01 --is_paired
上一步就已经生成结果了,这一步只是格式化,下面的代码会筛选出HLA-A基因的分型结果
perl parse_result.pl Allelelist.txt result.txt | grep "^A\*" | sort -k2 -n -r > HLA.txt
格式化之后的结果,内容如下
A*01:01:01:01 17.4022266628604A*11:01:01 12.0376819868684
共两列,第一列为Allel, 第二列为该Allel区域的平均测序深度。
联系客服