需求客户的一个简单需求:
示例文件gene.txt image.png test.vcf image.png 代码实现run.sh cat $1 |while read gene chr from todo#echo $chr $from $toif echo $2 |grep -q '.*.vcf.gz$';thenvcftools --gzvcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out $gene.$chr.$from-$to elif echo $2 |grep -q '.*.vcf$';thenvcftools --vcf $2 --chr $chr --from-bp $from --to-bp $to --recode --recode-INFO-all --out $gene.$chr.$from-$tofidone 运行 生成结果: image.png 补充说明以上代码中利用了vcftools工具,以及shell中读取每行文件的每个字段进行赋值。 vcftools还能提取某个具体位置的SNP: vcftools --gzvcf test.vcf.gz --positions specific_position.txt --recode --out specific_position.vcf specific_position.txt文件格式如下: 1 842013
1 891021
1 903426
1 949654
1 1018704 除了vcftools,bcftools和plink等工具也能实现类似的功能。 bcftools filter test.vcf.gz --regions 9:4700000-4800000 > out.vcf 但bcftools要求vcf必须是gz格式,如不是,则需要进行转化(直接用gzip不行): bcftools view test.vcf -Oz -o test.vcf.gzbcftools index test.vcf.gz 需要格外注意的是,vcf中的染色体名称要和提取文件中的染色体名保持一致,如Chr1或chr1或1。 或者: bcftools view -S keep.list test.vcf >sub_indv.vcf keep.list可以是“染色体+具体位置”两列,也可以是“染色体+起始+终止”三列: chr1 27639
chr1 60383
chr2 60469
chr3 60516
chr4 60534#或者chr1 1 1000
chr1 2000 4500 在plink中,可以指定特定的样本(keep)或SNP(extract)。 指定样本提取: plink --bfile file --noweb --keep sampleID.txt --recode --make-bed --out sample sampleID.txt第一列为提取的样本Family ID,第二列为Within-family ID(IID)。 指定位点提取: plink --bfile file --extract snp.txt --make-bed --out snp
snp.txt文件中一个SNP名称一行。 |
|
来自: 刘得光3p6n6zqq > 《SNP》