EBI提供HLA序列BLAST

基于我们有的HLA序列,可以和HLA序列的数据库比较,看与哪个HLA allele最相似。

HLA (human leukocyte antigen,人类白细胞抗原)是人类主要组织相容性复合体(major histocompatibility complex,MHC)的表达产物,根据HLA抗原结构、功能及组织分布的不同,分为I类,II类,III类分子,其中I类分子包括HLA-A,-B,-C系列抗原,广泛分布于各组织有核系统表面。

BLAST表示局部比对搜索工具,用来将新的序列与已有的数据库中的序列进行比较,可以发现区域的相似性,进而为功能和进化研究提供线索。
EBI(欧洲生物信息学中心)提供基于IPD-IMGT/HLA(IMGT国际免疫遗传学数据库)数据库的BLAST库。BLAST工具会搜索数据库中的HLA allele的核苷酸、蛋白质及相关对的序列。

HLA BLAST在线服务的链接如下:
https://www.ebi.ac.uk/Tools/services/web_ncbiblast/toolform.ebi?tool=ncbiblast&context=nucleotide&database=imgthla

1)选择数据库


选择IMGT下的IMGT/HLA (genomic)数据据,表示与HLA的基因组序列进行比对,如果需要和编码区序列进行比对,可以选择IMGT/HLA (cds)。

2)输入待比对的序列

BLAST支持简并碱基的比对,比如用R表示G或A (puRine),用A表示M表示A或C。

3)设置参数


选择blastn表示进行核苷酸序列比对,task可以用megablast也可以用blastn,其他参数可以选择默认。

4)提交任务

点击submit提交任务。

5)等待结果

一般情况下,会在10秒左右进入比对结果页面。结果页面包含相似性和显著性结果。

测试序列

CGAACTTGGCGGGTCTCAGCCCTCCTSGCCCCAGGCTCCCACTCCATGAGGTATTTCTACACCGCCATGTCCCGGCCCGG
CCGCGGGGAGCCCCGCTTCATCRCMGTGGGCTACGTGGACGACACSCAGTTCGTGMGGTTCGACAGCGACGCCRCGAGTC
CGAGRRKGGMGCCSCGGGCGCCRTGGRTRGAGCAGGAGGGGCCGGAGTATTGGGACCGGAASACACAGATCTMCAAGGCC
MASGCACAGACTGACCGAGAGAACCTGCGSAACCTGCTCGGCTACTACAACCAGAGCGAGGCCGGTGAGTGACCCCGGCC
CGGGGCGCAGGTCACGACTCCCCATCCCCCACGKACGGCCCGGGTCGCCCCGAGTCTCCGGGTCCGAGATCCGCCTCCCT
GAGGCCGCGGGACCCGCCCAGACCCTCGACCGGCGAGAGCCCCAGGCGCGTTTACCCGGTTTCATTTTCAGTTGAGGCCA
AAATCCCCGCGGGTTGGTCGGGGCGGGGCGGGGCTCGGGGGACGGGGCTGACCGCGGGGCCGGGGCCAGGGTCTCACACT
TGGCAGACGATGTATGGCTGCGACCTGGGGCCGGACGGGCGCCTCCTCCGCGGGCATAACCAGTTAGCTACGACGGCAGG

该序列的测试结果,与本地blast的结果一致。

参考:

https://www.ebi.ac.uk/ipd/imgt/hla/blast.html

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

将VCF文件中的突变拆分成SNP和INDEL

VCFTOOLS

得到SNP

vcftools --vcf X.vcf --remove-indels --out X.snps --recode --recode-INFO-all

得到INDEL

vcftools --vcf X.vcf --keep-only-indels --out X.indel --recode --recode-INFO-all

GATK

得到SNP

java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R reference.fasta \
    -V X.vcf   \
    -selectType SNP \
    -o X.snps.vcf

得到INDEL

java -jar GenomeAnalysisTK.jar \
    -T SelectVariants \
    -R reference.fasta \
    -V X.vcf  \
    -selectType INDEL \
    -o X.indel.vcf

本着有轮子不造轮子的原则,可以用VCFTOOLS和GATK来实现,当然如果想自己拆分的话,可以根据VCF中是否有SNP和INDEL的tag标签,或者根据ALT和REF中的碱基长度是否一致来实现拆分。

参考:
https://www.biostars.org/p/48204/
https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

分析带UMI标签的测序数据

分析带UMI标签的测序数据

检测癌组织的低频突变,为了提高检测低频突变的灵敏度,往往进行高深度的测序。但样本之间存在交叉污染,测序有存在一定概率的错误,这些因素会导致高深度测序过程中将假阳性的信号放到,得到假阳性的结果。解决交叉污染的方法,有公司比如IDT采用唯一配对的样本index,只有配对的index中的reads才属于特定样本。解决测序错误的方法,研究人员在建库的时候,先对分子加上UMI碱基,unique molecular identifier -> UMI,然后根据来源于同一个分子的测序数据进行测序错误修正,得到正确的分子序列。两种方法结合可以减少交叉污染提高准确性(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5759201/)。

如图中所示,左侧一个分子被测了5次,其中第二次有一个测序错误,但该错误并没有在每个测序数据中出现,所以在后续合成一个分子的时候,测序错误被修正,只保留了真正的突变。(https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5852328/)

常规的肿瘤配对测序分析,或者遗传性突变位点的分析,并不需要UMI信息,所以包含UMI的数据分析是需要不一样的分析流程来得到准确的分析结果,其中包括提取UMI分子标签,合并来自同一个分子的测序reads,低频突变检测而非胚系突变检测等。

大致流程为:

-----prepare analysis ready BAM file------
|         FASTQ
|            ↓
|          uBAM
|            ↓  Extract UMI
|          uBAM
|             ↓  Align uBAM and merge 
|           BAM
-----call consensus------
|             ↓  Group Reads By Umi
|           BAM
|             ↓  Group Reads By Umi
|           BAM
|             ↓  Call Molecular Consensus Reads
|           uBAM
|             ↓  Align uBAM and merge
|           BAM
|             ↓  Filter Consensus Reads
|           BAM
|             ↓  Clip
|           BAM
-----Vardict------
|             ↓  Call
|            VCF

一,得到包含UMI分子标签信息的BAM文件

UMI信息,应该从fastq中的配对的reads中提取,但fastq不能存储更多的信息,所以需要先将fastq转成uBAM文件,提取uBAM文件中的UMI分子标签信息,将该信息通过RX标签写入uBAM文件中,通过uBAM和BAM文件合并,把RX信息合并到比对到的BAM文件中进行下一步分析。

1)生成uBAM

java -Xmx8G -jar picard.jar FastqToSam \
    FASTQ=$fq1    \
    FASTQ2=$fq2   \
    OUTPUT=test.uBAM \
    READ_GROUP_NAME=test \
    SAMPLE_NAME=test \
    LIBRARY_NAME=test \
    PLATFORM_UNIT=HiseqX10  PLATFORM=illumina \
    RUN_DATE=`date --iso-8601=seconds`

2)提取UMI信息

java -jar fgbio.jar ExtractUmisFromBAM \
    --input=test.uBAM  \
    --output=test.umi.uBAM \
    --read-structure=2M148T 2M148T \
    --single-tag=RX \
    --molecular-index-tags=ZA ZB

此时,uBAM文件中RX标签记录着UMI的信息,2M148T表示前两个碱基是UMI分子标签,148个template 测序序列,配对read在uBAM文件中有如下信息

R1----ZA:Z:TC ZB:Z:TA RG:Z:test       RX:Z:TC-TA
R2----ZA:Z:TC ZB:Z:TA RG:Z:test       RX:Z:TC-TA

继续阅读

什么是ubam文件,为什么ubam文件比fastq文件好

ubam是Unmapped BAM Format,是BAM文件的一个变种,里面的read是未经map的。大部分测序供应FASTQ文件,这是最常见的测序分析的起始文件。FASTQ文件的优势是,压缩比比bam文件好,解压速度快。

但与ubam文件相比,FASTQ并不是最理想的:

1)单个文件更容易分析

在双端测序中,有些软件希望配对的reads放在一个文件中,有些希望配对的文件,有些软件直接根据read在文件中的位置判断read是否配对,当然现在FASTQ通过在文件中加入/1和/2来表明配对的read解决这一问题。在生信分析的时候,一个FASTQ文件往往和另一个FASTQ文件关联配对,比如R1和R2,往往会花费更多时间来验证read是否配对。但如果通过bam文件的话,会更简单。只需要在FLAG这个地方加入对应的值,比如77和141,就能指定。单个文件更加简单,能储存更多metadata的信息。

2)BAM文件可以储存更多信息

实际上,正是这个原因,我才开始用ubam文件的。
因为FASTQ文件格式的问题,很多其他信息不能存在FASTQ文件中,而BAM/SAM文件可以储存很多信息,只需加tag或者header就行。如果你想放,可以放用了什么软件,平台,样本标签等。
上图是不利用ubam文件的流程,下图是利用ubam文件的流程。

这个我体会很深,比如我在分析带分子标签UMI的数据时,需要在bam文件中添加read的UMI标签序列,具有相同的UMI标签的序列说明来自同一个分子。如果我将带UMI标签的reads直接比对基因组,UMI标签序列会影响比对质量,而且后续对bam文件的处理提取UMI时,需要考虑CIGAR标签,这个时候就头大了。实际上,我是先将fastq转成ubam文件(picard的FastqToSam工具),然后提取ubam文件中的UMI,生成新的ubam,将ubam转成fastq文件(reads都放在一个文件fastq文件中,bwa能识别其中的配对reads),与参考基因组进行比对。这个时候问题就来了,比对之后生成的bam文件,不包含我提取的UMI标签序列信息,因为这个信息在ubam文件中。不要着急,picard的另外一个工具MergeBamAlignment,可以将ubam文件和bam文件合并,这样ubam文件中的tag信息,比如UMI标签信息,就会补充到新的bam文件中。这样生成的bam文件供下一步数据分析。

3)BAM能实现随机读写

BAM文件采用的是BGZF压缩,可以通过随机读写快速读取文件特定位置的数据。对于未比对的reads(unaligned reads),这个功能不是很吸引人,因为后续序列比对都是顺序比对的,不可能跳着比对。

参考:
https://blastedbio.blogspot.com/2011/10/fastq-must-die-long-live-sambam.html
https://gatkforums.broadinstitute.org/gatk/discussion/5990/what-is-ubam-and-why-is-it-better-than-fastq-for-storing-unmapped-sequence-data
https://gatkforums.broadinstitute.org/gatk/discussion/11694/why-is-converting-from-fastq-to-ubam-nesessary-before-preprocessing
https://gatkforums.broadinstitute.org/gatk/discussion/6484/how-to-generate-an-unmapped-bam-from-fastq-or-aligned-bam

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################

fastq压缩之后的gzip文件大小与样本数据量

在测序的时候,我们先拿到的是样本fastq压缩后的gzip文件,这个时候可能最关心的是数据量够不够,那么fastq.gz文件大小和测序数据量有什么关系呢。
我用Miseq测序数据(gz文件200M左右),Hiseq panel(gz文件50M左右)和WES测序数据(gz文件4G左右)进行了简单的分析。有意思的地方是,虽然R1和R2的数据量是一样的,解压出来的文件大小是一样的,但R2的gzip文件总比R1大。不管是Miseq还是Hiseq的panel测序,压缩后的R2均大于R1文件,且文件越小,差异越大。

因为R1和R2数据量相同的原因,我只看R1的真实文件和gz文件大小与数据量之间的关系。
数据量=FASTQ文件行数/4*151/1000/1000 单位为M
真实文件大小估计=FASTQ文件行数/4*357/1024/1024 单位为M,预测值,差别不大,因为FASTQ文件中每四行357个字符(和平台和设置有关系),每个字符1byte。
GZ文件大小通过ll -h查看
因为FASTQ文件是规范的,每四行字符基本一致,所以FASTQ真实文件大小和数据量成正比。比如我前面提到的每四行有357个字符,其中序列只占151个字符,也就是说FASTQ文件大小大概是测序量的357/151≈2.3倍多。但因为FASTQ文件为文本文件,占用空间较大,所以一般将FASTQ文件压缩成gzip格式文件。

用gzip文件大小除以估计的真实文件大小,得到压缩比,发现压缩比和测序平台有关系,Miseq的压缩比在0.35左右,Hiseq平台中的panel测序在0.24左右,WES在0.2左右。可以发现WES的数据的压缩效率最高,猜测不同平台在压缩过程中设置的压缩比例不同。Hiseq测序平台的gzip文件压缩性能大于Miseq平台。

以Miseq的gzip文件大小和数据量作图,可以很直观的发现gzip文件大小和数据量之间的线性关系。如图,斜率是1.245,截距可以忽略不计,因为纵坐标很大。也就是说R1文件的数据量约等于gzip文件的1.245倍。

从另外一个角度计算,R1文件的数据量=gzip文件大小/压缩效率/2.3,则因为Miseq的压缩效率在0.35,Hiseq在0.23左右,所以Miseq R1文件的数据量大概是gzip文件的1.242倍(1/0.35/2.3),和上图拟合的结果很相近,Hiseq约在1.89(1/0.23/2.3)倍。真实的数据量其实是R1文件数据量的2倍(R2文件中的数据量和R1文件相同)。

结论:不同平台和不同测序方法之间的fastq文件压缩比例并不一致,但同一种方法和平台内的压缩比例是相近的,因而可以根据gzip文件的大小推测出测序的数据量。
我手头的数据中,Miseq的测序量约是R1 gzip文件大小的2.49倍,Hiseq的测序量约是R1 gzip文件大小的3.78倍。不同平台和方法之间,这个值是不固定的,可以根据手头的数据计算一下这个比例,就能迅速的估算出样本的数据量。

#####################################################################
#版权所有 转载请告知 版权归作者所有 如有侵权 一经发现 必将追究其法律责任
#Author: Jason
#####################################################################