01从vcf中挑选指定区间或指定样本

一.指定区间

1.0项目应用

1.list文件没有头信息 直接处理
$ vcftools --vcf samples.pop.snp.anno.result.list --chr Chr4A --from-bp 742765746 --to-bp 744656845  --recode --recode-INFO-all --out Chr4A.742765746-744656845.snp.anno.result.list # snp
$ vcftools --vcf samples.pop.indel.anno.result.list --chr Chr4A --from-bp 742765746 --to-bp 744656845  --recode --recode-INFO-all --out Chr4A.742765746-744656845.indel.anno.result.list # indel

2.vcf有head信息 直接处理,头信息会包含所有染色体
$ vcftools --gzvcf/--vcf samples.pop.snp.anno.result.vcf.gz/samples.pop.snp.anno.result.vcf --chr Chr4A --from-bp 742765746 --to-bp 744656845  --recode --recode-INFO-all --out Chr4A.742765746-744656845.snp.anno.result.vcf 

3.对于vcf文件 需要进一步处理 头信息只包含所需染色体
$ mv all.sv.gt.anno.result.vcf unsorted.all.sv.gt.anno.result.vcf # sv
$ bcftools sort unsorted.all.sv.gt.anno.result.vcf -o all.sv.gt.anno.result.vcf # sv

$ cp samples.pop.snp.anno.result.vcf samples.pop.snp.anno.result.vcf.bak # snp
$ cp samples.pop.indel.anno.result.vcf samples.pop.indel.anno.result.vcf.bak # indel
$ cp all.sv.gt.anno.result.vcf all.sv.gt.anno.result.vcf.bak # sv

$ bgzip samples.pop.snp.anno.result.vcf # 压缩 snp
$ bgzip samples.pop.indel.anno.result.vcf # 压缩 indel
$ bgzip all.sv.gt.anno.result.vcf # 压缩 sv

$ bcftools index samples.pop.snp.anno.result.vcf.gz # 索引 snp
$ bcftools index samples.pop.indel.anno.result.vcf.gz # 索引 indel
$ bcftools index all.sv.gt.anno.result.vcf.gz # 索引 sv

$ (bcftools view -h samples.pop.snp.anno.result.vcf.gz | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h samples.pop.snp.anno.result.vcf.gz | grep -w "ID=Chr4A")" ; bcftools view -H  -r Chr4A:742765746-744656845   samples.pop.snp.anno.result.vcf.gz  ) | bcftools view -Oz -o Chr4A.742765746-744656845.snp.anno.result.vcf.gz # 输出为.gz文件

$ (bcftools view -h samples.pop.snp.anno.result.vcf.gz | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h samples.pop.snp.anno.result.vcf.gz | grep -w "ID=Chr4A")" ; bcftools view -H  -r Chr4A:742765746-744656845   samples.pop.snp.anno.result.vcf.gz  ) > Chr4A.742765746-744656845.snp.anno.result.vcf # 直接输出 snp
$ (bcftools view -h samples.pop.indel.anno.result.vcf.gz | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h samples.pop.indel.anno.result.vcf.gz | grep -w "ID=Chr4A")" ; bcftools view -H  -r Chr4A:742765746-744656845   samples.pop.indel.anno.result.vcf.gz  ) > Chr4A.742765746-744656845.indel.anno.result.vcf # 直接输出 indel
$ (bcftools view -h all.sv.gt.anno.result.vcf.gz | grep -v "##contig" | sed "/^##reference/a $(bcftools view -h all.sv.gt.anno.result.vcf.gz | grep -w "ID=Chr4A")" ; bcftools view -H  -r Chr4A:742765746-744656845   all.sv.gt.anno.result.vcf.gz  ) > Chr4A.742765746-744656845.sv.gt.anno.result.vcf # 直接输出 sv

$ gunzip -c samples.pop.snp.anno.result.vcf.gz > samples.pop.snp.anno.result.vcf # 保留源文件 解压 snp
$ gunzip -c samples.pop.indel.anno.result.vcf.gz > samples.pop.indel.anno.result.vcf # 保留源文件 解压 indel
$ gunzip -c all.sv.gt.anno.result.vcf.gz > all.sv.gt.anno.result.vcf # 保留源文件 解压 sv

1.1 vcftools

区间位置文件

chr_position.txt（基因染色体区间坐标文件）

chr1    15645891        15646070
chr1    15646256        15646447
chr1    15646534        15646644
chr1    15646772        15646811
chr1    15646861        15646899
chr1    15647179        15647279
chr1    15647377        15647555

test.vcf（snp信息vcf文件，不一定需要vcf文件，只要满足下面示例文件第一二列坐标就可以提取）

#CHROM  POS        ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  CF2
chr01    112       94      .       T       C       33689.1 .       AC=1;AF=0.472;AN=2;BaseQRankSum=-0.076;DP=3953;ExcessHe
chr02    1123665   547     .       A       G       67292.8 .       AC=1;AF=0.489;AN=2;BaseQRankSum=0;DP=4273;ExcessHet=129
chr05    112553    1245    .       A       T       69849.7 .       AC=1;AF=0.466;AN=2;BaseQRankSum=-0.339;DP=4231;ExcessHe
chr05    142554    2086    .       C       G       80488.5 .       AC=1;AF=0.506;AN=2;Base
chr07    1111125   2334    .       A       G       42490.1 .       AC=1;AF=0.466;AN=2;BaseQRankSum=-0.028;DP=3806;ExcessHe
chr10    1111445   2335    .       A       T       53023.1 .       AC=1;AF=0.466;AN=2;BaseQRankSum=0.199;DP=3835;ExcessHet

脚本如下

cat $1 | while read chr from to
do
    echo $chr $from $to
    # Extract base filename without extension from $2
    base=$(basename $2)
    prefix="${base%.*}"
    
    if echo $2 | grep -q '.*.vcf.gz
```bash
bash Chr_position_tk_snp.sh          位置文件          输入文件

参数说明

vcftools --gzvcf input.vcf --chr n --recode – recode-INFO-all --stdout | gzip -c > output.vcf.gz  
说明：  
–gzvcf：处理压缩格式的vcf文件（可替换为–vcf）  
–chr n：选择染色体n，例：–chr 1  
–recode：重新编码为vcf文件,有过滤操作都要加上--recode  
–recode-INFO-all：将输出的文件保存所有INFO信息  
–stdout：标准输出，后接管道命令  
–gzip -c：压缩

提取某个具体位置的SNP

vcftools --vcf test.vcf --positions specific_position.txt --recode --out specific_position.vcf

specific_position.txt文件格式如下：

1.2 bcftools需要tbi索引（vcf必须是gz格式）

提取位置或区间

1.直接输入
bcftools filter test.vcf.gz --regions 9:4700000-4800000 > out.vcf

bcftools view -r NC_061100.1:10224668-10230294 variants_vqsr_filter/samples.gatk.055.con.snp.vqsr.filter.vcf.gz|less -S # 验证ok

2.大批量取位点
染色体及位置文件chr_pos.list,示例如下
chr1    27639
chr1    60383
chr2    60469
chr3    60516
chr4    60534  
或者直接给区间也行  
chr1　　1　　1000  
chr1　　2000　　4500

bcftools view -R chr_pos.list samples.snp.region.filter.vcf.gz > new.samples.snp.region.filter.vcf

压缩建索引

bcftools要求vcf必须是gz格式，如不是，则需要进行转化（直接用gzip不行）

## 方法1
bgzip samples.snp.region.filter.vcf && tabix -p vcf samples.snp.region.filter.vcf.gz

## 方法2
bcftools view test.vcf -Oz -o test.vcf.gz
bcftools index test.vcf.gz

ref：[bcftools 提取vcf（snp/indel）文件子集 - 天使不设防 - 博客园]
(https://www.cnblogs.com/mmtinfo/p/11945592.html)

1.3 plink

可以指定特定的样本（keep)或SNP(extract)。

指定位点提取

plink --bfile file --extract snp.txt --make-bed --out snp

snp.txt文件中一个SNP名称一行。

1.4 参考

Ref：https://www.cnblogs.com/chenwenyan/p/9151672.html
https://blog.csdn.net/weixin_34387468/article/details/94519445
https://www.cnblogs.com/mmtinfo/p/11945592.html
https://www.cnblogs.com/chenwenyan/p/8991417.html
如何从vcf文件中批量提取一系列基因的SNP位点？

二.指定样品

2.1 方法1: 使用vcftools

vcftools --vcf input.vcf --keep samples.txt --recode --out output

samples.txt文件应包含每行一个样本名。

2.2 方法2: 使用bcftools

bcftools view -s sample1,sample2,sample3 input.vcf > output.vcf

或者使用样本文件(每行一个样本名):

bcftools view -S samples.txt input.vcf > output.vcf

2.3 方法3: 使用GATK SelectVariants

gatk SelectVariants \
  -V input.vcf \
  -O output.vcf \
  -sn sample1 -sn sample2 -sn sample3

2.4 plink

plink --bfile file --noweb --keep sampleID.txt --recode --make-bed --out sample

sampleID.txt第一列为提取的样本Family ID，第二列为Within-family ID（IID）。

; then
vcftools --gzvcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out ${prefix}_$ {chr}. $f r o m -$ {to}
elif echo $2 | grep -q '.*.vcf

参数说明

提取某个具体位置的SNP

{{CODE_BLOCK_6}}
specific_position.txt文件格式如下：

1.2 bcftools需要tbi索引（vcf必须是gz格式）

提取位置或区间

压缩建索引

bcftools要求vcf必须是gz格式，如不是，则需要进行转化（直接用gzip不行）
{{CODE_BLOCK_9}}
ref：[bcftools 提取vcf（snp/indel）文件子集 - 天使不设防 - 博客园]
(https://www.cnblogs.com/mmtinfo/p/11945592.html)

1.3 plink

可以指定特定的样本（keep)或SNP(extract)。

指定位点提取

{{CODE_BLOCK_10}}
snp.txt文件中一个SNP名称一行。

1.4 参考

Ref：https://www.cnblogs.com/chenwenyan/p/9151672.html
https://blog.csdn.net/weixin_34387468/article/details/94519445
https://www.cnblogs.com/mmtinfo/p/11945592.html
https://www.cnblogs.com/chenwenyan/p/8991417.html
如何从vcf文件中批量提取一系列基因的SNP位点？

二.指定样品

2.1 方法1: 使用vcftools

{{CODE_BLOCK_11}}
samples.txt文件应包含每行一个样本名。

2.2 方法2: 使用bcftools

{{CODE_BLOCK_12}}
或者使用样本文件(每行一个样本名):

2.3 方法3: 使用GATK SelectVariants

2.4 plink

{{CODE_BLOCK_15}}
sampleID.txt第一列为提取的样本Family ID，第二列为Within-family ID（IID）。

; then
vcftools --vcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out ${prefix}_$ {chr}. $f r o m -$ {to}
else
vcftools --vcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out ${prefix}_$ {chr}. $f r o m -$ {to}
fi
done

{{CODE_BLOCK_4}}

>[!note] 参数说明
{{CODE_BLOCK_5}}
### 提取某个具体位置的SNP
{{CODE_BLOCK_6}}
specific_position.txt文件格式如下：
{{CODE_BLOCK_7}}
## 1.2 bcftools需要tbi索引（vcf必须是gz格式）
### 提取位置或区间
{{CODE_BLOCK_8}}
### 压缩建索引
bcftools要求vcf必须是gz格式，如不是，则需要进行转化（直接用gzip不行）
{{CODE_BLOCK_9}}
ref：[bcftools 提取vcf（snp/indel）文件子集 - 天使不设防 - 博客园]
(https://www.cnblogs.com/mmtinfo/p/11945592.html)
## 1.3 plink
可以指定特定的样本（keep)或SNP(extract)。
### 指定位点提取
{{CODE_BLOCK_10}}
snp.txt文件中一个SNP名称一行。
## 1.4 参考
> Ref：[https://www.cnblogs.com/chenwenyan/p/9151672.html](https://www.cnblogs.com/chenwenyan/p/9151672.html)  
> [https://blog.csdn.net/weixin_34387468/article/details/94519445](https://blog.csdn.net/weixin_34387468/article/details/94519445)  
> [https://www.cnblogs.com/mmtinfo/p/11945592.html](https://www.cnblogs.com/mmtinfo/p/11945592.html)  
> [https://www.cnblogs.com/chenwenyan/p/8991417.html](https://www.cnblogs.com/chenwenyan/p/8991417.html)
[如何从vcf文件中批量提取一系列基因的SNP位点？](https://www.cnblogs.com/jessepeng/p/14530891.html)

# 二.指定样品
## 2.1 方法1: 使用vcftools
{{CODE_BLOCK_11}}
`samples.txt`文件应包含每行一个样本名。
## 2.2 方法2: 使用bcftools
{{CODE_BLOCK_12}}
或者使用样本文件(每行一个样本名):
{{CODE_BLOCK_13}}
## 2.3 方法3: 使用GATK SelectVariants
{{CODE_BLOCK_14}}
## 2.4 plink
{{CODE_BLOCK_15}}
sampleID.txt第一列为提取的样本Family ID，第二列为Within-family ID（IID）。