1.Linux服務(wù)器安裝Aspera

Note:如果服務(wù)器無(wú)法用其中一個(gè)版本的Aspera，建議更換Aspera的版本加叁。
參考:
http://www.reibang.com/p/ff826d6591f5 --服務(wù)器1參考(相對(duì)可以用較高版本)
https://github.com/indexofire/bac-ngs-book/blob/master/chapter_1/sra.md --服務(wù)器2參考(相對(duì)只能用低版本)

2.SNP數(shù)據(jù)下載——PRJEB32865為例(EBI上下載)

#Note:先從ENA上下載aspera地址(見(jiàn)下圖)----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)分個(gè)體。(3)已經(jīng)合并好。
#case 1,2:bcftools合并
# 6.output file name:這里我以物種名和project兩級(jí)目錄存放數(shù)據(jù)党饮。
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] 部分下載下來(lái)可能沒(méi)有index
ls *.vcf.gz > VCF_LIST
for j in $(cat VCF_LIST)
do
        echo "*** ${i} begin to build index ***"
        #bcftools index -t $i
        tabix -p vcf $i 
        echo "*** Finished ***"
done
## 7.2 merge
### 7.2.1 分染色體的VCF文件合并
bcftools concat --allow-overlaps --threads 10 -f VCF_LIST -Oz > $File_Name_Input
### 7.2.2 分個(gè)體的VCF文件合并
bcftools merge *.vcf.gz -Oz --force-samples -o $File_Name_Input

分染色體

分個(gè)體

3.vcf文件轉(zhuǎn)bed,bim,fam

# 8.vcf--bed,bim,fam
plink2 --vcf  ${File_Name_Input} --max-alleles 2 --chr-set 90 --make-bed --allow-extra-chr --out ${File_Name_Output}
#最大的染色體編號(hào)(可以通過(guò)下面方式獲取常染色體條數(shù))
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)報(bào)錯(cuò)1:Invalid chromosome code ’chr29‘
人類有23對(duì)染色體，PLINK默認(rèn)把23解釋為x染色體,因此對(duì)于非人類物種(超過(guò)22條常染色體)笆焰，需要用--chr-set指定常染色體數(shù)目劫谅，這里不妨把--chr-set設(shè)大點(diǎn)。
(2)報(bào)錯(cuò)2:xx.bim cannot contain multiallelic variants
設(shè)定--max-alleles 2嚷掠。
(3)報(bào)錯(cuò)3:vcf file decompression failure: Malformed BGZF block
vcf文件解壓失敗(少部分會(huì)出現(xiàn)這種情況)捏检，主要是因?yàn)橄螺d下來(lái)文件不全，通過(guò)在腳本中添加md5check解決不皆。
(4)報(bào)錯(cuò)4:Invalid chromosome code ’JH373233.1‘
設(shè)定--allow-extra-chr
(5)報(bào)錯(cuò)5:Failed to open xx.vcf.gz: could not load index
說(shuō)明下載下來(lái)的vcf文件沒(méi)有對(duì)應(yīng)的index贯城，需要用bcftools構(gòu)建index。
(6)報(bào)錯(cuò)6:Duplicate sample names (xxx), use --force-samples to proceed anyway
bcftools合并過(guò)程中報(bào)錯(cuò)霹娄，加上--force-samples參數(shù)即可能犯。
例:A.vcf.gz含有樣本S1、S2犬耻、S3踩晶，B.vcf.gz含有樣本S3,S4，使用force-samples將會(huì)出現(xiàn)4個(gè)樣本S1,S2,2;S3,S4枕磁。
(7)染色體編號(hào)問(wèn)題:多倍體以及部分SNP數(shù)據(jù)中染色體編號(hào)可能不是正常的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文件。
(8)報(bào)錯(cuò)7:Error: Line xx of --vcf file has a GT half-call.
解決方法設(shè)置--vcf-half-call 'missing' 

報(bào)錯(cuò)1

報(bào)錯(cuò)2

報(bào)錯(cuò)3

報(bào)錯(cuò)4

報(bào)錯(cuò)5

報(bào)錯(cuò)6

報(bào)錯(cuò)7