SNPsplit—区分等位基因reads神器

SNPsplit

  区分等位基因reads神器  

基因组印记（Genomic imprinting）一直是二倍体哺乳动物的研究热点，特别是在印记基因这一块，但往往由于大家对这方面的软件了解较少，因此总是心有余而力不足，不知如何下手！今天，小编就给大家介绍一款十分简单、好使的软件—SNPsplit。
说到基因组印记，就不得不提到两个专有名词，一个是等位基因特异表达(allele specific expression)、另一个是等位基因特异甲基化(allele specific methylation)。那么如何从全基因组水平上，较为全面定义子代继承亲本的印记模式呢？这是一个需要解决的问题。
2016年，Babraham生物信息组织发布了一款专门用来区分亲本来源reads的软件SNPsplit,它通过SAM/BAM文件reads上覆盖的已知SNP位点信息，能够将reads分配给其中一个等位基因。

1. 软件说明

SNPsplit软件只需要提供用来区分印记来源的SNP信息（通过双亲本VCF文件得到，详见下述实操部分——4.1、4.2），就可以针对生信常用软件（包括Bowtie2, TopHat, STAR, HISAT2, HiCUP 和Bismark）比对后的bam区分其亲本来源。

该软件的主要包括以下2部分模块：1）SNPsplit_genome_preparation重新建立比对N-masked基因组，其原理如下图，举例流程详见下述4.2部分。

2）SNPsplit
针对重新比对的bam文件，区分亲本来源reads。举例流程详见4.3、4.4。

2. 软件安装

SNPsplit是用Perl编写的，在github下载包后，解压就可以直接使用了。注意SNPsplit软件需要利用samtools工具。github上该软件的所在路径如下：https://github.com/FelixKrueger/SNPsplit

tar -zxvf SNPsplit-0.3.4.tar.gz
cd SNPsplit-0.3.4
#测试是否能运行
./SNPsplit_genome_preparation --help
./SNPsplit --help

我们主要是使用以下两个命令：

1

SNPsplit_genome_preparation

2

SNPsplit

3. 候选SNP位点的筛选

按照实操4.2，便可生成区分F1子代的候选SNP位点文件：

all_MotherM_SNPs_FatherL_reference.based_on_GRCm38.txt（双亲本纯合且不同）

当然为了得到更加可信的结果，我们可以对这些SNP位点进一步过滤，比如筛选reads深度都>=10的SNP位点。此部分可以根据自己项目研究的实际情况来做。

4. 软件实操

那么接下来，小编以转录本数据为例，给大家分享下利用SNPsplit软件定义印记基因的流程。

4.1 数据准备

双亲本merged SNP位点的VCF文件：merged.vcf基因组原始参考序列：genome.faF1子代转录组原始数据：test_R1.fq.gz，test_R2.fq.gz

#使用GATK合并双亲本SNP位点vcf文件
java -XX:ParallelGCThreads=4 -Xmx18g -jar GenomeAnalysisTK.jar -T CombineVariants -R genome.fa -V FatherL.snp.vcf.gz  -V MotherM.snp.vcf.gz  -o merged.vcf 

（代码稍长，动动手指往右滑动）

4.2 生成新的参考基因组

利用SNPsplit_genome_preparation生成所需要的N-masker基因组和候选SNP位点。

mkdir genome
mv genome.fa genome
SNPsplit_genome_preparation --vcf_file merged.vcf --strain FatherL  --dual_hybrid --strain2 MotherM --reference_genome genome --genome_build GRCm38

（代码稍长，动动手指往右滑动）
SNPsplit_genome_preparation会生成以下结果文件：

其中红色框中的部分是接下来区分genome1（父本）,genome2（母本）所需要的依赖文件。

cd FatherL_MotherM_dual_hybrid.based_on_GRCm38_N-masked
cat *fa > all_N-masked.GRCm38.N-masked.fa

N-masked.fa：all_N-masked.GRCm38.N-masked.fa
候选SNP位点：all_MotherM_SNPs_FatherL_reference.based_on_GRCm38.txt

4.3 重新比对

首先，我们需要用上一步骤得到的N-masked.fa参考基因组根据相应的比对软件重新生成所需的新的index。之后按照以下命令重新比对。注意SNPsplit需要按照read ID排序的bam文件作为输入。

hisat2 -x /HISAT2_new_index/all_N-masked.GRCm38.N-masked --rf --rna-strandness RF --known-splicesite-infile splicesites.txt --dta -t -p 4 -1 test_R1.fq.gz -2 test_R2.fq.gz  --un test --no-softclip| samtools test/test.all.bam -S -

samtools sort -n -@ 4 -m 2G test.all.bam -o test.sortn.bam

（代码稍长，动动手指往右滑动）

4.4 区分genome1/2 reads

SNPsplit --snp_file all_MotherM_SNPs_FatherL_reference.based_on_GRCm38.txt --paired --no_sort -o test --singletons test.sortn.bam

（代码稍长，动动手指往右滑动）

我们会得到以下结果文件：

test.sortn.genome1.bam #父本来源reads
test.sortn.genome2.bam #母本来源reads

4.5 定义印记基因

最后，我们只需要将genome1和genome2来源的reads当成2个不同条件来源的比较组（多个样本可以当成replicates），用htseq-count统计所有注释基因的count数，然后用DEseq2做差异分析便可。
由于印记基因是只在表达双亲本一方的等位基因，因此这里我们可以选择比较高的cutoff，小编这里推荐，qval<0.05, abs(log2FC)>=3。
按照以上流程定义出印记基因后，我们还可以与数据库中的结果进行比对。数据库路径如下：http://www.geneimprint.com/site/genes-by-species

5. 参考

1.Krueger F, Andrews SR. SNPsplit: Allele-specific splitting of alignments between genomes with known SNP genotypes. F1000Res. 2016 Jun 23;5:1479. doi: 10.12688/f1000research.9037.1http://europepmc.org/article/MED/274297432.https://github.com/FelixKrueger/SNPsplit3.https://github.com/FelixKrueger/SNPsplit/blob/master/SNPsplit_User_Guide.md

作者 | Jenny
审稿 | 童蒙
编辑 | angelica

我们有信心与武汉共渡难关关注阿拉丁，为科研实用而生

往期回顾

姿势已摆好

就等你点啦