引言
基因组印记（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
我们主要是使用以下两个命令：
- SNPsplit_genome_preparation
- 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.fa
F1子代转录组原始数据：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. 参考
- 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.1
http://europepmc.org/article/MED/27429743 - https://github.com/FelixKrueger/SNPsplit
- https://github.com/FelixKrueger/SNPsplit/blob/master/SNPsplit_User_Guide.md
作者 | Jenny
审稿 | 童蒙
编辑 | angelica
我们有信心与武汉共渡难关!!! 关注阿拉丁，为科研实用而生!
往期回顾
转录组数据定量归一化python3字体解决大挖掘
SV碰上三维结构，探究Duplication通过影响TAD如何引发疾病？
神灯宝典之三代重测序分析实录（二）
神灯宝典之PB三代重测序分析实录（一）
从生到死，是什么驱动了染色质的分相变化？
生信老司机教你做基因组项目
我命由天不由我！ecDNA,你所不知道的癌症内幕
如何使用软件自动对变异进行ACMG打分
网友评论