師兄讓下載bcftools膀息,學(xué)習(xí)一下勉痴。
1.下載
電腦太不給力了,wget都沒有盔粹,用brew下載一個
太慢了隘梨,在github直接下載。
https://github.com/samtools/bcftools/releases
安裝看readme
共享目錄出錯舷嗡,
這個時候就要在/etc/ld.so.conf中加入xxx.so所在的目錄轴猎。
一般而言,有很多so檔會在/usr/local/lib這個目錄下
在/etc/ld.so.conf中加入/usr/local/lib這一行后,執(zhí)行 在etc中執(zhí)行l(wèi)dconfig 命令來更新生效进萄。
./bftools 使用
2.samtools語句
介紹samtools和bcftools的語句與用法 :http://www.reibang.com/p/da32e3c3a168
samtools基本語句和功能:
1.view 查看bam捻脖,sam文件的內(nèi)容
2.sort用來對bam文件進行排序
3.merge將2個或2個以上的已經(jīng)sort的bam文件合并成一個bam文件锐峭,合并后的文件不需要再次sort,是已經(jīng)sort過了的可婶。
4.index對bam文件建立索引沿癞,生成后綴為.bai的文件,用于快速檢索reads矛渴。需要對bam文件先進行排序椎扬,否則會報錯。很多時候都需要索引文件的存在具温,特別是顯示序列比對的情況下蚕涤。比如samtools tview,gbrowse2等。
5.faidx對fasta序列建立索引铣猩,生成后綴為.fai的文件揖铜。該命令也能依據(jù)索引文件快速提取fasta文件中的某一條序列。
6.tview能直觀的顯示出reads比對到基因組的情況达皿,和基因組瀏覽器有點類似天吓。
7.flagstat給出bam文件的比對結(jié)果的summary。
8.depth得到每個位點的測序深度峦椰,并輸出到標準輸出
samtools depth a.bam | less
9.rmdup去除PCR重復(fù)和光學(xué)重復(fù)龄寞,重復(fù)的reads僅保留一對。
10.reheader替換bam文件的頭文件
11.cat鏈接多個bam文件们何,適用于非sorted的bam文件
12.mpileup(現(xiàn)在也屬于bcftools)
使用samtools的子命令mpileup分析參考序列上的每個堿基位點的比對結(jié)果萄焦,并生成VCF/BCF格式文件,BCF是VCF的二進制文件冤竹。再使用Bcftools對VCF/BCF格式文件進行SNP/Indel calling拂封。其中,Bcftools是附屬于samtools的程序鹦蠕。
后面接Bcftools進行variation calling
bcftools進行SNP calling:? www.reibang.com/p/af93ca26ba44
time $bcf mpileup ${precallbam} --fasta-ref ${reference} | $bcf call -mv -o ./out/bcfcall.vcf && echo " ** bcf snp calling done **"
在進行SNP calling 時冒签,必須選擇一種算法,有兩種calling算法可供選擇钟病,分別對應(yīng)-c和-m參數(shù)萧恕。-c參數(shù)對應(yīng)consensus-caller算法,-m參數(shù)對應(yīng)multiallelic-caller算法肠阱,后者更適合多種allel和罕見變異的calling票唆。
-v參數(shù)也是常用參數(shù),作用是只輸出變異位點的信息屹徘,如果一個位點不是snp/indel, 不會輸出走趋。
3.bam文件格式
samtools view -H SRR1770413.sorted.markdup.bam
詳細講解bam文件:http://www.reibang.com/p/364e640d3c9f
4.bftools合并vcf文件
師兄用bftools合并vcf文件,這里的vcf文件不是gatk得到的噪伊,沒有combine過簿煌。
http://www.reibang.com/p/59728f20d38c
5 vcf 文件格式
http://www.reibang.com/p/38f734ae47f5
http://www.reibang.com/p/957efb50108f
剛剛看到一個實戰(zhàn):
http://www.reibang.com/p/6f3198b7a070
發(fā)現(xiàn)gatk在github上有workflow: https://github.com/gatk-workflows
其中vcf文件生成步驟:
1.HaplotyperCaller? ? ? 2.MergeVcfs? ? ?3.GenomicsDBImport? ? ?4.GenotypeGVCFs
在github上選擇workflow打開氮唯,查看wdl文件中的command了解流程