1.Linux服务器安装Aspera
Note:如果服务器无法用其中一个版本的Aspera，建议更换Aspera的版本。
参考:
https://www.jianshu.com/p/ff826d6591f5 --服务器1参考(相对可以用较高版本)
https://github.com/indexofire/bac-ngs-book/blob/master/chapter_1/sra.md --服务器2参考(相对只能用低版本)
2.SNP数据下载——PRJEB32865为例(EBI上下载)
#Note:先从ENA上下载aspera地址(见下图)----filereport_analysis_PRJEB32865_tsv.txt
## 1.remove the first col+separate the second and third col
awk '{$1="";print $0}' file*.txt | awk '{print $2}' > aspera.txt
awk '{$1="";print $0}' file*.txt | awk '{print $1}' > md5.txt
## 2.remove first row+replace ; with \n
sed -i '1d' aspera.txt
sed -i '1d' md5.txt
sed -i "s/;/\n/g" aspera.txt
sed -i "s/;/\n/g" md5.txt
## 3.merge aspera+md5
awk -F "/" '{print $NF}' aspera.txt > file_name.txt
paste md5.txt file_name.txt > md5_check.txt
## 4.remove repeat(for download)
sort -n aspera.txt | uniq > aspera_uniq.txt
sort -n md5_check.txt | uniq > md5_check_uniq.txt
## 5.Ensure full download
for i in $(cat aspera_uniq.txt)
do
echo "*** ${i} begin to download ***"
ascp -v -k 1 -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@${i} ./
File_Name=$(echo $i | awk -F "/" '{print $NF}')
grep -w "${File_Name}$" md5_check_uniq.txt > tmp.md5
md5sum -c tmp.md5 > tmp.md5check
# MD5 check pass
while true
do
#Continue,until download entirely
grep "FAILED" tmp.md5check > /dev/null
if [ $? -eq 0 ]
then
#delete+downlaod again
rm $File_Name
echo "*** ${i} begin to download ***"
ascp -v -k 1 -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@${i} ./
File_Name=$(echo $i | awk -F "/" '{print $NF}')
grep -w "${File_Name}$" md5_check_uniq.txt > tmp.md5
md5sum -c tmp.md5 > tmp.md5check
else
echo "*** Finished ***"
break
fi
done
done
ENA下载aspera地址
#Note:下载得到vcf文件可能有三种情况:(1)分染色体。(2)分个体。(3)已经合并好。
#case 1,2:bcftools合并
# 6.output file name:这里我以物种名和project两级目录存放数据。
PWD=$(pwd)
Species=$(echo $PWD | awk -F '/' '{print $(NF-1)}')
Project=$(echo $PWD | awk -F '/' '{print $NF}')
File_Name_Input=$Species"_"$Project".vcf.gz"
File_Name_Output=$Species"_"$Project
# 7.merge(optional)
# 7.1 Index(tbi)[-f:overwrite index,-t:*.tbi] 部分下载下来可能没有index
ls *.vcf.gz > VCF_LIST
for j in $(cat VCF_LIST)
do
echo "*** ${i} begin to build index ***"
bcftools index -t $i
echo "*** Finished ***"
done
## 7.2 merge
bcftools merge *.vcf.gz -Oz --force-samples -o $File_Name_Input
分染色体
分个体
3.vcf文件转bed,bim,fam
# 8.vcf--bed,bim,fam
plink2 --vcf ${File_Name_Output} --max-alleles 2 --chr-set 90 --make-bed --allow-extra-chr --out ${File_Name_Output}
#最大的染色体编号(可以通过下面方式获取常染色体条数)
Max_Chr=$(cat ${File_Name_Output}.bim | awk '{print $1}' | awk ' !x[$0]++' | grep '^[0-9]' | awk 'NR==1{max=$1;next}{max=max>$1?max:$1}END{print max}')
echo "*** ${Species} has ${Max_Chr} autosomal chromosomes ***"
Note:
(1)报错1:Invalid chromosome code ’chr29‘
人类有23对染色体，PLINK默认把23解释为x染色体,因此对于非人类物种(超过22条常染色体)，需要用--chr-set指定常染色体数目，这里不妨把--chr-set设大点。
(2)报错2:xx.bim cannot contain multiallelic variants
设定--max-alleles 2。
(3)报错3:vcf file decompression failure: Malformed BGZF block
vcf文件解压失败(少部分会出现这种情况)，主要是因为下载下来文件不全，通过在脚本中添加md5check解决。
(4)报错4:Invalid chromosome code ’JH373233.1‘
设定--allow-extra-chr
(5)报错5:Failed to open xx.vcf.gz: could not load index
说明下载下来的vcf文件没有对应的index，需要用bcftools构建index。
(6)报错6:Duplicate sample names (xxx), use --force-samples to proceed anyway
bcftools合并过程中报错，加上--force-samples参数即可。
例:A.vcf.gz含有样本S1、S2、S3，B.vcf.gz含有样本S3,S4，使用force-samples将会出现4个样本S1,S2,2;S3,S4。
(7)染色体编号问题:多倍体以及部分SNP数据中染色体编号可能不是正常的chr1,chr2,chr3/1,2,3----需要自己用bcftools修改染色体名。
bcftools annotate --rename-chrs chr_name_conv.txt *.vcf.gz -Oz -o $File_Name_Input
#Note:
chr_name_conv.txt:第一列为待命名染色体，第二列为最终想命名成的染色体名称。
*.vcf.gz:待重命名的vcf文件。
$File_Name_Input:输出vcf文件。
报错1
报错2
报错3
报错4
报错5
报错6
网友评论