如何從vcf檔案中批量提取一系列基因的SNP位點？

目錄客戶的乙個簡單需求：

我有一批功能基因位點，想從重測序的群體材料中找到這些位點，如何批量快速獲得？

執行sh run.sh gene.txt test.vcf，或sh run.sh gene.txt test.vcf.gz生成結果：

以上**中利用了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.gz

bcftools 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名稱一行。

ref：

利用Python批量重新命名一系列檔名雜亂的檔案

假設目錄下面有這樣一系列命令雜亂的檔案 openfoam training part 1.pdf openfoam training part 2.pdf openfoam training part 3 pdf 不僅序號被放在最後，而且還有許多多餘的空格。現在批量將這些檔案重新命名，去掉 並把序號...

如何複製陣列中一系列元素的元素

本例項主要介紹如何使用 array 類的copy 方法來複製陣列中一系列的元素。copy 方法從指定的源索引開始，複製 array 中的一系列元素，將他們貼上到另乙個 array 中 從指定的目標索引開始 長度和索引指定為 64位整數。其方法有多種過載形式，本例項所使用的過載形式如下 public ...

R語言讀取Excel檔案的一系列陷阱

你想用r讀取乙個excel檔案，你覺得這事沒啥難的，就像所有的檔案讀取，只需要知道檔名就萬事大吉了。於是，你把1.xls放到讀取.r的資料夾下面，重新命名為1.csv，開啟rstudio，執行下面這條語句 a出現了下面的報錯 error in file file,rt cannot open the...