首先在linux 里配置conda
下載
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive/Anaconda3-5.3.1-Linux-x86_64.sh
chmod +x Anaconda3-5.3.1-Linux-x86_64.sh
bash Anaconda3-5.3.1-Linux-x86_64.sh
安裝完畢,如果忘記選擇yes,敲conda命令報錯“command not found" 加上source /root/anaconda3/etc/profile.d/conda.sh
conda env list 得到 /root/anaconda3
export PATH=~/anaconda3/bin:$PATH
不然全局無法使用conda命令,(但是重啟putty好像就不管用了,還不清楚原因)
vs code可以不安裝
安裝tree命令,yum install tree
tree -af 可以查看樹形文件結構
snakemake是純python的任務流程工具(基於python3),以前商業環境用過control-M
https://snakemake.readthedocs.io/en/stable/
首先做個變異檢測,就是和標准的序列做對比,有點類似於代碼的compare,用過Beyond Compare或svn和git的筒子們應該很熟悉了,
但是基因序列是個非常大的序列文件,在linux也沒有windows那樣簡便的圖形操作見面,而且,這種對比工作是大量重復的,需要腳本化。
cd snakemake-snakemake-tutorial-623791d7ec6d
conda env create --name snakemake-tutorial --file environment.yaml
--------------------------------------------------------------
export PATH=~/anaconda3/bin:$PATH
source activate snakemake-tutorial
--------------------------------------------------------------
# 退出當前環境
source deactivate
這里使用到Samtools工具,具體使用方法可以參考https://blog.csdn.net/g863402758/article/details/53081342
他是一個用於處理sam與bam格式的工具軟件,能夠實現二進制查看、格式轉換、排序及合並等功能,
結合sam格式中的flag、tag等信息,還可以完成比對結果的統計匯總。同時利用linux中的grep、awk等操作命令,
還可以大大擴展samtools的使用范圍與功能。
conda install snakemake
conda install samtools
bowtie2和samtools都是對比工具,bowtie2暫時沒安裝,安裝方法先記錄下
sudo wget https://jaist.dl.sourceforge.net/project/bowtie-bio/bowtie2/2.3.4.1/bowtie2-2.3.4.1-linux-x86_64.zip
unzip bowtie2-2.3.4.1-linux-x86_64.zip
vi /etc/environment
添加 bin 目錄的路徑,並用 :
隔開
source /etc/enviroment 使配置生效
開始寫job腳本
rule bwa_map:
input:
"data/genome.fa", "data/samples/A.fastq" output: "mapped_reads/A.bam" shell: """ bwa mem {input} | samtools view -Sb - > {output} """
期間一直出一個錯誤,說Command must be given as string after the shell keyword
運行snakemake -np mapped_reads/A.bam
檢查一下是否會出錯
執行這個job,把-n去掉
可以看到,生成了A.bam文件
rule bwa_map:
input:
"data/genome.fa", "data/samples/{sample}.fastq" output: "mapped_reads/{sample}.bam" shell: """ bwa mem {input} | samtools view -Sb - > {output} """
將A改成{sample},在輸入命令的時候加上你的參數,自動匹配上了,(注意此時文件夾貌似只能有一個腳本文件),cp了一個好像報錯了
接下來,要做排序了,代碼最后一起貼
可以使用dag
選項和dot
命令對“規則的執行和依賴關系”進行可視化,
snakemake --dag sorted_reads/{A,B,C}.bam.bai | dot -Tpdf > dag.pdf 這個命令好像會報錯
snakemake --dag sorted_reads/{A,B,C}.bam.bai | dot -Tsvg > dag.svg
整合之前的BAM文件,做基因組變異識別
SAMPLES=["A","B","C"] rule bcftools_call: input: fa="data/genome.fa", bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES), bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES) output: "calls/all.vcf" shell: "samtools mpileup -g -f {input.fa} {input.bam} | " "bcftools call -mv - > {output}"
其中expand是自動匹配變量求文件路徑的語法糖
檢查一下,snakemake -np calls/all.vcf

最后出report,以上都是在規則里執行shell腳本,snakemake的一個優點就是可以在規則里面寫Python腳本,只需要把shell改成run,此外還不需要用到引號。
測試一下,snakemake -np report.html
畫出流程圖
snakemake --dag report.html | dot -Tsvg > final.svg
執行一下:snakemake -p report.html
可以看到生成了報告文件

到此,還有
rule all:
log:
多線程thread:
-j 指定cpu核心
params:
加載configfile: "config.yaml"
這幾個功能沒有操作,留個以后有空再處理
最后,在新建一個snakemake項目時,都先用conda create -n 項目名 python=版本號創建一個全局環境,用於安裝一些常用的軟件,例如bwa、samtools、seqkit等。然后用如下命令將環境導出成yaml文件
conda env export -n 項目名 -f environment.yaml
以后再部署的時候,
只需要conda env create -f environment.yaml
這個過程類似於ghost系統,或者打包虛擬機類似
參考了以下網址,感謝!
https://www.jianshu.com/p/8e57fd2b81b2
http://pedagogix-tagc.univ-mrs.fr/courses/ABD/practical/snakemake/snake_intro.html