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

最新推荐文章于 2024-06-13 07:03:11 发布

子鹿学生信

最新推荐文章于 2024-06-13 07:03:11 发布

阅读量1.3k

点赞数

文章标签：程序人生

本文链接：https://blog.csdn.net/weixin_42192188/article/details/134084795

版权

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

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

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

这是我的第36篇文章。

以下是正文：

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

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

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

1. 下载小鼠参考数据集

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

网页链接见UCSC：https://hgdownload.soe.ucsc.edu/downloads.html#mouse

 wget --timestamping 
        'ftp://hgdownload.cse.ucsc.edu/goldenPath/mm10/chromosomes/*'
# 可以解压合并到一起
gunzip <file>.fa.gz
cat *fa > GRCm38.fa

# GRCm39参考基因组
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/635/GCA_000001635.9_GRCm39

下载的小鼠参考基因组在用于比对前，需要先进行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}/GRCm38.fa.gz -O GRCm38.fa.dict

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

wget --recursive --no-parent --no-directories \
--accept vcf*vcf.gz \
ftp://ftp.ncbi.nih.gov/snp/organisms/archive/mouse_10090/VCF/

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

# SNP
# NCIB的资源
wget ftp://ftp.ncbi.nih.gov/snp/organisms/archive/mouse_10090/VCF/00-All.vcf.gz 
wget ftp://ftp.ncbi.nih.gov/snp/organisms/archive/mouse_10090/VCF/00-All.vcf.gz.tbi

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

# Sanger Mouse Genetics Programme (Sanger MGP)的资源
wget ftp://ftp-mouse.sanger.ac.uk/REL-1505-SNPs_Indels/mgp.v5.merged.snps_all.dbSNP142.vcf.gz

# Indel
wget ftp://ftp-mouse.sanger.ac.uk/REL-1505-SNPs_Indels/mgp.v5.merged.indels.dbSNP142.normed.vcf.gz \
-O mgp.v5.indels.vcf.gz

# illumina也提供了很多资源
http://igenomes.illumina.com.s3-website-us-east-1.amazonaws.com/Mus_musculus/UCSC/mm10/Mus_musculus_UCSC_mm10.tar.gz

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

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

# keep only passing and append 
zcat mgp.v5.indels.vcf.gz | 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

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

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

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

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

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

# 删除非染色体信息的vcf：vcf_chr_NotOn.vcf等
rm vcf_chr_NotOn.vcf.gz vcf_chr_Multi.vcf.gz vcf_chr_Un.vcf.gz vcf_chr_AltOnly.vcf.gz  -rf

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

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

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

http://ftp.cbi.pku.edu.cn/pub/CGPS_download/
https://ftp.ncbi.nih.gov/snp/organisms/archive/mouse_10090/VCF/genotype/SC_MOUSE_GENOMES.genotype.vcf.gz
http://crispor.tefor.net/genomes/mm10/
https://www.sanger.ac.uk/data/mouse-genomes-project/
https://ftp.ebi.ac.uk/pub/databases/mousegenomes/REL-1505-SNPs_Indels/
http://ftp.ensembl.org/pub/release-106/variation/gvf/mus_musculus/
https://support.illumina.com/sequencing/sequencing_software/igenome.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}_1.fq.gz -I $rawdata_path/${f}_2.fq.gz -o $output_path/${f}.R1.fq.gz -O $output_path/${f}.R2.fq.gz
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/GRCm38.fa.gz

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 "@RG\tID:QSY${name}\tLB:QSY${name}\tSM:${name}\tPL:ILLUMINA" \
       ${ref} ${rawdata_path}/${name}.R1.fq.gz ${rawdata_path}/${name}.R2.fq.gz \
       | 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}.align.reorder.bam \
       SEQUENCE_DICTIONARY=${ref_dict}

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

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

$GATK MarkDuplicates \
 -I ${output_path}/${name}.align.reorder.sorted.bam \
 -M ${output_path}/${name}.markdup_metrics.txt \
 -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/C57BL_6NJ.mgp.v5.indels.dbSNP142.normed.vcf \
 --known-sites $GATK_bundle/C57BL_6NJ.mgp.v5.snps.dbSNP142.vcf \
 -O $output_path/${sample}.sorted.markdup.recal_data.table && echo "** ${sample}.sorted.markdup.recal_data.table done **" 

time $GATK ApplyBQSR \
 --bqsr-recal-file $output_path/${sample}.sorted.markdup.recal_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}.HC.vcf.gz && echo echo "** ${sample}.HC.vcf.gz 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}.HC.vcf.gz \
 --select-type-to-include SNP \
 -O $Variants_path/${sample}.SNP.vcf.gz

$gatk SelectVariants \
 -R $ref \
 -V $Variants_path/${sample}.HC.vcf.gz \
 --select-type-to-include INDEL \
 -O $Variants_path/${sample}.INDEL.vcf.gz
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--

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

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

本文由 mdnice 多平台发布

子鹿学生信

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
3
评论
必收藏！超详细的小鼠全基因组测序数据分析流程！

下载的vcf文件中header中的##reference=GCF_000001635.24，一定要存在该reference，就是说你用到的reference 文件名字一定要是header中的这个名字。我目前在构建一个能够互相分享系统深入学生信的社群（免费提供内部生信学习资料），且大家也可以通过创作获得收益，有认同我理念的小伙伴可以来找我聊聊哈！如果是人类样本，一般用GATK的bundle资源进行VQSR，但是小鼠的资源非常少，而且样本数少也不建议用VQSR，因此后面只能自己设置参数进行硬质控过滤了！
复制链接

扫一扫