必收藏！超详细的小鼠全基因组测序数据分析流程！

嗨！你好，我是超速成长的子鹿，致力于构建一个欢乐好学有深度且能分享运营收益的生信社群​！

在INFJ和ENFP中反复横跳，不给自己设限的六边形梦想实干家。

更多信息请到文章末尾查看！​

这是我的第36篇文章。

以下是正文：

不知道为什么，网上很少有小鼠全基因组测序数据分析流程，可能是因为小鼠的基因组资源比较少吧。

这里我汇总了目前经常用到的小鼠参考数据库资源，辅助进行变异检测。

同时给了一个可行的小鼠全基因组分析流程供大家参考。

1. 下载小鼠参考数据集

1. 首先下载小鼠的 GRCm38参考基因组

网页链接见UCSC：.html#mouse

 wget --timestamping 
        '/*'
# 可以解压合并到一起
gunzip <file>.fa.gz
cat *fa > GRCm38.fa

# GRCm39参考基因组
wget .9_GRCm39

2. 下载的小鼠参考基因组在用于比对前，需要先进行index和dict处理：

samtools=~/miniconda/envs/WGS/bin/samtools
ref=/home/GRCm38.fa
$samtools faidx $ref

GATK=~/software/gatk-4.3.0.0/gatk
Java=/usr/bin/java
GRCm38_path=/home/GRCm38

${GATK} CreateSequenceDictionary -R ${GRCm38_path}/ -O GRCm38.fa.dict

3. 下载dbSNP数据库的 GRCm38 VCF文件（这个文件是按染色体分开的）

wget --recursive --no-parent --no-directories 
--accept  
/

后面又找到合并好的文件，包括SNP和Indel：

# SNP
# NCIB的资源
wget . 
wget .bi

其实Sanger Mouse Genetics Programme (Sanger MGP)和illumina也提供了很多资源

# Sanger Mouse Genetics Programme (Sanger MGP)的资源
wget .v5.merged.snps_all.

# Indel
wget .v5.merged. 
-O mgp.v5.

# illumina也提供了很多资源
.

后面发现上面的参考文件是全部的位点，需要进一步过滤出PASS的位点

# take header first
zcat mgp.v5. | head -1000 | grep "^#" | cut -f 1-8 
> mgp.v5.indels.pass.chr.vcf

# keep only passing and append 
zcat mgp.v5. | grep -v "^#" | cut -f 1-8 
| grep -w "PASS"  >> mgp.v5.indels.pass.chr.vcf

# 排序
gatk SortVcf -SD GRCm38_68.dict -I mgp.v5.indels.pass.chr.vcf -O mgp.v5.indels.pass.chr.sort.vcf

# rm .idx
# rm mgp.v5.indels.pass.chr.sort.vcf.idx

4. 给染色体数字前面加上"chr"，这步是保证跟参考基因组的染色体编号一样

回去看了下之前比对用的参考基因组，确实是带chr的！

# 修改vcf文件，加上chr
for vcf in $(ls -1 *.) ; do
  vcf_new=${vcf/./.vcf}
  echo $vcf
  zcat $vcf | sed 's/^([0-9XY])/chr1/' > $vcf_new
  rm -fv $vcf
done

# 修改头文件中的染色体编号
sed -i 's/##contig=<ID=/##contig=<ID=chr/g' d.dbSNP142.head1000.vcf

# 文件较大的话用这个代码，用sed只修改前1000行
sed '1,1000s/##contig=<ID=/##contig=<ID=chr/g' d.dbSNP142.vcf > d.dbSNP142.sub.vcf
# 也可以用AWK
awk 'NR<=1000 {gsub("##contig=<ID=","##contig=<ID=chr");} 1' d.dbSNP142.vcf > d.dbSNP142.vcf

# 删除非染色体信息的vcf：vcf_chr_NotOn.vcf等
rm vcf_chr_ vcf_chr_ vcf_chr_ vcf_chr_  -rf

*需要注意的是： 下载的vcf文件中header中的##reference=GCF_000001635.24，一定要存在该reference，就是说你用到的reference 文件名字一定要是header中的这个名字。GATK对header的要求比较严格。

下载的vcf文件中染色体的开始要和基因组参考序列fasta文件的一致。*

其他小鼠相关资源链接如下：

1. /
2. .
3. /
4. /
5. /
6. /
7. .html

2. 原始数据过滤

#!/bin/bash

rawdata_path=/dssg/home/Rawdata
output_path=/dssg/home/01FqQC

for f in KO-299-LFJ9356_L3 WT-298-LFJ9355_L3; do 
fastp -i $rawdata_path/${f}_ -I $rawdata_path/${f}_ -o $output_path/${f}. -O $output_path/${f}.
done

3. bwa比对到参考基因组

rawdata_path=/home/01FqQC
output_path=/home/02Mapping

bwa=~/miniconda/envs/WGS/bin/bwa
samtools=~/miniconda/envs/WGS/bin/samtools
bamtools=~/miniconda/envs/WGS/bin/bamtools
picard=~/miniconda/envs/WGS/bin/picard
GATK=~/software/gatk-4.3.0.0/gatk
Java=/usr/bin/java

ref=/home/

for name in KO-299-LFJ9356_L3; do 

echo -e The Step MappingAndMarkDup started at `date` "n" >>${output_path}/run.log;

bwa mem -t 20 -M -R "@RGtID:QSY${name}tLB:QSY${name}tSM:${name}tPL:ILLUMINA" 
       ${ref} ${rawdata_path}/${name}. ${rawdata_path}/${name}. 
       | gzip -3 > ${output_path}/${name}.align.sam

time $samtools view -@ 20 -bS ${output_path}/${name}.align.sam 
        -o ${output_path}/${name}.align.bam

${picard} ReorderSam 
       INPUT=${output_path}/${name}.align.bam 
       OUTPUT=${output_path}/${name}.der.bam 
       SEQUENCE_DICTIONARY=${ref_dict}

 ${samtools} sort -@ 20 -m 8G 
        ${output_path}/${name}.align.bam 
        -o ${output_path}/${name}.der.sorted.bam

 ${samtools} index -@ 20 ${output_path}/${name}.der.sorted.bam

$GATK MarkDuplicates 
 -I ${output_path}/${name}.der.sorted.bam 
 -M ${output_path}/${name}. 
 -O ${output_path}/${name}.sorted.markdup.bam && echo "** ${name}.sorted.bam MarkDuplicates done **"

time $samtools index ${output_path}/${name}.sorted.markdup.bam && echo "** ${name}.sorted.markdup.bam index done **"

time $GATK BaseRecalibrator 
 -R $ref 
 -I $output_path/${sample}.sorted.markdup.bam 
 --known-sites $GATK_bundle/p.v5.d.vcf 
 --known-sites $GATK_bundle/p.v5.snps.dbSNP142.vcf 
 -O $output_path/${sample}.al_data.table && echo "** ${sample}.al_data.table done **" 

time $GATK ApplyBQSR 
 --bqsr-recal-file $output_path/${sample}.al_data.table 
 -R $ref 
 -I $output_path/${sample}.sorted.markdup.bam 
 -O $output_path/${sample}.sorted.markdup.BQSR.bam && echo "** ApplyBQSR done **"

## 为${sample}.sorted.markdup.BQSR.bam构建Index，这是继续后续步骤所必须的
time $samtools index $output_path/${sample}.sorted.markdup.BQSR.bam && echo "** ${sample}.sorted.markdup.BQSR.bam index done **"

4. 生成vcf文件

# 生成vcf文件
ref=/home/GRCm38.fa
GATK_bundle=/home/mm10_GATK/Indel

for sample in KO-299-LFJ9356_L3 WT_298_LFJ9355_L3; do
time $GATK HaplotypeCaller 
        -R $ref 
        -I $output_path/${sample}.sorted.markdup.BQSR.bam 
        -O $Variants_path/${sample}. && echo echo "** ${sample}. done ** "
done

# 将SNP和INDEL分开
for sample in KO-299-LFJ9356_L3 WT_298_LFJ9355_L3; do
# sample=KO-299-LFJ9356_L3

$gatk SelectVariants 
 -R $ref 
 -V $Variants_path/${sample}. 
 --select-type-to-include SNP 
 -O $Variants_path/${sample}.

$gatk SelectVariants 
 -R $ref 
 -V $Variants_path/${sample}. 
 --select-type-to-include INDEL 
 -O $Variants_path/${sample}.
done

5. Filtering Variants

分别对WT和KO样本的SNP和Indel进行硬质控过滤

如果是人类样本，一般用GATK的bundle资源进行VQSR，但是小鼠的资源非常少，而且样本数少也不建议用VQSR，因此后面只能自己设置参数进行硬质控过滤了！

# Filtering SNPs

rawdata_path=/home/01FqQC
output_path=/home/02Mapping
Variants_path=/home/03Variants

bwa=~/miniconda/envs/WGS/bin/bwa
samtools=~/miniconda/envs/WGS/bin/samtools
bamtools=~/miniconda/envs/WGS/bin/bamtools
picard=~/miniconda/envs/WGS/bin/picard
gatk=~/software/gatk-4.3.0.0/gatk
Java=/usr/bin/java

ref=/home/GRCm38.fa

for sample in KO-299-LFJ9356_L3 WT_298_LFJ9355_L3; do
# sample=KO-299-LFJ9356_L3

$gatk VariantFiltration 
 -R $ref 
 -V $Variants_path/${sample}.SNP.vcf 
 --filter-expression "QD < 2.0" --filter-name "QD" 
 --filter-expression "SOR > 3.0" --filter-name "SOR" 
 --filter-expression "FS > 60.0" --filter-name "FS" 
 --filter-expression "MQ < 40.0" --filter-name "MQ" 
 --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum" 
 --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum" 
 -O $Variants_path/${sample}.SNP.filtered.vcf
done


for sample in KO-299-LFJ9356_L3 WT_298_LFJ9355_L3; do
# sample=KO-299-LFJ9356_L3

$gatk VariantFiltration 
 -R $ref 
 -V $Variants_path/${sample}.INDEL.vcf 
 --filter-expression "QD < 2.0" --filter-name "QD" 
 --filter-expression "SOR > 3.0" --filter-name "SOR" 
 --filter-expression "FS > 60.0" --filter-name "FS" 
 --filter-expression "MQ < 40.0" --filter-name "MQ" 
 --filter-expression "MQRankSum < -12.5" --filter-name "MQRankSum" 
 --filter-expression "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum" 
 -O $Variants_path/${sample}.INDEL.filtered.vcf
done

统计各个条件下位点数量，WT通过PASS为328130个，KO为1166109个

grep -v "^#" WT_298_LFJ9355_L3.INDEL.filtered.vcf | cut -f 7 | sort | uniq -c

grep -v "^#" KO-299-LFJ9356_L3.SNP.filtered.vcf | cut -f 7 | sort | uniq -c

过滤出VCF文件中PASS的位点

awk -F 't' '{if($0 ~ /#/) print; else if($7 == "PASS") print}' WT_298_LFJ9355_L3.SNP.filtered.vcf > WT_298_LFJ9355_L3.SNP.PASS.vcf

awk -F 't' '{if($0 ~ /#/) print; else if($7 == "PASS") print}' KO-299-LFJ9356_L3.INDEL.filtered.vcf > KO-299-LFJ9356_L3.INDEL.PASS.vcf

合并SNP和Indel文件

gatk=~/software/gatk-4.3.0.0/gatk
Variants_path=/home/03Variants

for sample in KO-299-LFJ9356_L3 WT_298_LFJ9355_L3; do
$gatk MergeVcfs 
 --INPUT $Variants_path/${sample}.SNP.PASS.vcf 
 --INPUT $Variants_path/${sample}.INDEL.PASS.vcf 
 --OUTPUT $Variants_path/${sample}.PASS.vcf
done

尝试直接比较两个VCF位点差异

# Output file comparing the sites in two vcf files
vcftools=~/miniconda/envs/WGS/bin/vcftools
$vcftools --vcf KO-299-LFJ9356_L3.PASS.vcf 
--diff WT_298_LFJ9355_L3.PASS.vcf 
--not-chr chr5_JH584297_random 
--diff-site --out KO_v_WT 
--diff-site-discordance --out KO_v_WT

7. 提取vcf中目标位置变异

vcftools=~/miniconda/envs/WGS/bin/vcftools
#$vcftools --vcf KO-299-LFJ9356_L3.PASS.vcf --chr chr7 
#  --out KO.Srcap.chr7
$vcftools --vcf WT_298_LFJ9355_L3.PASS.vcf --chr chr7 --from-bp 127513018 --to-bp 127561219 --recode --out wt_Sr
$vcftools --vcf KO-299-LFJ9356_L3.PASS.vcf --chr chr7 --from-bp 127513018 --to-bp 127561219 --recode --out KO_Sr

慢慢理解世界，慢慢更新自己。共勉～

--THE END--

我目前在构建一个能够互相分享系统深入学生信的社群（免费提供内部生信学习资料），且大家也可以通过创作获得收益，有认同我理念的小伙伴可以来找我聊聊哈！更多系统学生信的教程可以直接加入社群获取！

欢迎大家私聊我，和子鹿做朋友~