单行命令
2021-11-03 本文已影响0人
琴酒martini
利用plink进行数据预处理（修剪SNP集）
！！！注意：如果要同时使用，正确做法是先SNP（geno）后样本(mind)
SNP缺失率多少用不了 - CSDN
1.缺失率 --geno 0.05
删除基因型缺失率大于5%的SNPs，大于95%的个体都具有的变异位点才保留
就是缺失，没测出来
2.最小等位基因频率 --maf 0.01
将次要等位基因频率maf < 0.01的SNP筛选出来并过滤掉，仅包括MAF >= 0.01的SNP，默认值为0.01（防止假阳性）
（1.2在近交系数、NJ树、iqtree、iBS、iBD、ML树前用）
3.LD过滤 --indep-pairwise 50 25 0.2
窗宽50、删除LD大于0.2的SNP对中一个、每次将窗口向前移动25个SNP
（3只在PCA 、ADMIXTURE前用）
芯片数据不需要ld过滤！！！！！
plink过滤.png
根据id提取vcf（提取重测序数据不会有日志等啰嗦行）
nohup /home/software/bcftools/bcftools view -S id.txt test.vcf.gz -Ov > sampleid.vcf &
# id.txt为个体名
根据染色体提取vcf
nohup /home/software/vcftools/src/cpp/vcftools --gzvcf test.vcf.gz --recode --recode-INFO-all --chr n --out out.vcf &
根据位置提取vcf（pos文件含染色体号和位置）
nohup /home/software/vcftools/src/cpp/vcftools --gzvcf test.vcf.gz --positions ./position.txt --recode --out test.vcf.gz-position &
根据区间提取vcf
nohup /home/software/bcftools/bcftools view test.vcf -Oz -o test.vcf.gz &
nohup /home/software/bcftools/bcftools index test.vcf.gz &
nohup /home/software/bcftools/bcftools filter test.vcf.gz --regions 18:14700001-14800001 > 14700001-14800001.vcf &
提取vcf文件中的个体list （不会影响fid、iid）
nohup /home/software/bcftools/bcftools query -l test.vcf > test-id.txt &
格式转化（不会影响fid、iid）
### vcf to map（此步骤只适合染色体号为123格式，如为NC，则map文件第一列为未知，即0）
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --plink --out test &
### map to bed
nohup /home/software/plink --allow-extra-chr --chr-set 29 --file test --make-bed --out test &
### bed to vcf
nohup /home/software/plink --allow-extra-chr --chr-set 29 --bfile test --recode vcf-iid --out test &
### vcf to map (重测序数据使用)
nohup /home/software/plink --vcf test.vcf --recode --allow-extra-chr --chr-set 30 --out test &
### map to vcf (用于beagle)
nohup plink --allow-extra-chr --chr-set 30 --file test --recode vcf-iid --out test &
去除多等位基因，indel
nohup /home/software/bcftools/bcftools view -m 2 -M 2 --type "snps" test.vcf.gz -O z -o test.record.snps.vcf.gz &
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --remove-indels (--max-missing 0.9 --maf 0.03) --min-alleles 2 --max-alleles 2 --recode --recode-INFO-all --out test.miss.snp.recode.vcf &
去除0号染色体
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --recode --recode-INFO-all --not-chr 0 --out test_noinclude0 &
过滤掉indel，只保留snp:
nohup /home/software/vcftools/src/cpp/vcftools --remove-indels --recode --recode-INFO-all --vcf test.vcf --stdout > test.snp.vcf &
过滤掉snp，只保留indel：
nohup /home/software/vcftools/src/cpp/vcftools --keep-only-indels --recode --recode-INFO-all --vcf test.vcf --stdout > test.indel.vcf &
统计indel碱基数量：
bcftools stat xxx.vcf>xxx.vcf.stat
在vcf文件中删除一个个体
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --recode --recode-INFO-all --stdout --remove-indv XXX > out.vcf &
缺失率统计
## 按照位点统计
nohup /home/software/vcftools/src/cpp/vcftools --gzvcf test.vcf.gz --missing-site --out test.SNP_missing &
## 按照个体统计
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --missing-indv --out test.SNP_missing &
解压缩保留源文件
nohup gunzip -c test.vcf.gz > test.vcf &
统计vcf压缩文件共有行
gzip -dc test.vcf.gz|wc -l
循环更改脚本(例如将
2,$6)
sed -i 's/ $0/$2,$6/g' test.sh (-i 为在文件中修改)
sort排序与去重复
# 以第2列来排序，n为将文本作为数值解释
sort -n -k 2 test.bed > test.sorted.bed
# 去重复，去除文件中相同的行
sort -n -k 2 -u test.bed > test-u.bed
# -t $'/t'指定tab为分隔符，n为将第一列按照数值排列（1，2，3），g为按照科学计数法排列（1，2，3）
sort -t $'/t' -k 1n,1 -k 3g,3 test.bed > test.sorted.bed
sort -k 1,1d -k 2,2n 第一列按照字符排序（1，10，11，12），第二列按照数字排序）
sort更改空格与tab键
sed -i 's/\t/ /g' test.bim #去除vcf文件中的\t键
sed -i 's/ /\t/g' test.bim
去除以#号开头的注释行
sed -i '/^#/d' id.txt
去除文件前36行
sed -i '1,36d' id.txt
去除重复
#重复行应为相邻状态
uniq test.txt > test-u.txt
vi打开vcf文件后/寻找的
/find
excel
=VLOOKUP(A1,excelA1:D1,4,0) 取对应值
=LEFT(A1,2) 取A1的左侧2列
=RIGHT(A1,2) 取A1的右侧2列
=ROUNDDOWN（A1,0）去除小数点
=LOOKUP(9^9,B2:E2) 取行最后一列数据
=LEN（A1）统计单元格字数
=IF(MOD(A1,3)=0,A1,"") 删选3倍数
=IF(COUNTIF(J$1:J1,J1)=1,COUNTIF(J:J,J1),"") 统计某列某单元格出现得次数
去除多余染色体
nohup /home/software/vcftools/src/cpp/vcftools --vcf test.vcf --recode --recode-INFO-all --not-chr 0 --out test_noinclude0 &
去除空白行
sed '/^$/d' test.txt > test-kongbai.txt
excel表格中将1月1日替换为1/1
1.将所有数据改为文本格式
2.点击替换
替换.png
3.44197 的格式选择文本
1|1的格式不设定
4.全部替换即可
linux寻找共同点
sort 1.txt > 1.sort.txt
sort 2.txt > 2.sort.txt
comm -12 1.sort.txt 2.sort.txt > out.txt
提取出文件2中与文件1中的相同行
grep -f 1.txt 2.txt
根据第一列提取对应染色体
awk '{if($1==18) print $0}' test.txt > chr18.test.txt
芯片数据与重测序合并删除列
sed -i 's/ /\t/g' banbenposition.txt tab键
sed -e "s/^M//" banbenposition.txt > banbenposition_huiche.txt 去除回车符
cat Data_position.recode.vcf| awk -F " " '{$1=null;print}' > test_1.vcf
cat test_1.vcf| awk -F " " '{$1=null;print}' > test_12.vcf
paste -d "" banbenposition_huiche.txt test_12.vcf > Data_position.recode_banben.vcf
将已经替换染色体的前两列v1.0版本与去除前两列信息的v3.1版本vcf进行合并。
sed -i 's/ /\t/g' Data_position.recode_banben.vcf
寻找重合的区间
bedtools | 快速筛选重合区间 - 知乎 (zhihu.com)
nohup /home/software/bedtools2/bin/bedtools intersect -a a.txt -b b.txt -wao > ab.window.ihs.intersect.bed &
# 打开 a.txt
1 2419435 2419436
1 2419463 2419464
1 2419517 2419518
1 2419597 2419598
1 2419765 2419766
1 2420278 2420279
1 2420283 2420284
1 2420335 2420336
# 打开 b.txt
1 2100001 2200001
1 2300001 2400001
1 2400001 2500001
1 2500001 2600001
1 2600001 2700001
1 2700001 2800001
1 5200001 5300001
1 6600001 6700001
1 7600001 7700001
# 打开 ab.window.ihs.intersect.bed
1 2798804 2798805 1 2700001 2800001 1
1 2799204 2799205 1 2700001 2800001 1
1 2799911 2799912 1 2700001 2800001 1
1 2800978 2800979 . -1 -1 0
1 2801066 2801067 . -1 -1 0
1 2801680 2801681 . -1 -1 0
1 2801738 2801739 . -1 -1 0
1 2801784 2801785 . -1 -1 0
# 打开结果文件，前四列代表文件一里的区间，第5至8列代表文件一与文件二重合的区间，第九列代表重合的长度
# 文件一与文件二的重合长度为1
# .代表文件一中的区间d在文件二中未找到重合区间
