NextFlow学习记录--序列比对和评估

在分析NGS数据（fastq格式）之前，需要有参考基因组，bed文件，输入数据

这里构建一个流程，分为两步：

1. 将fastq比对为sam，并转成bam（需安装bwa或bowtie，还需安装samtools）
   

2. 对bam的某些区域进行测序深度评估，并绘图（需安装gnuplot）

首先用人类hg19的 Exons的bed文件，随机选取某些基因的Exons，作为测试bed。如果用整个bed，运算时间较长。

其次构建 nextflow的script文件，以nf文件为后缀

最后，nextflow run tutorial.nf 可以看到运行结果

如果尚未安装nextflow ， 采用以下命令安装

curl -s https://get.nextflow.io | bash

随机产生bed文件，有两种方法(这里选取了15行)

python shuf_bed.py 15 | awk '{print $1"\t"$2"\t"$3}'> test1.bed
## 或者
perl -MList::Util=shuffle -e 'print shuffle(<STDIN>);' < All_EXONhg19.bed | head -15 > test1.bed

shuf_bed.py 的内容如下：

import random
import sys
bedLines = open('All_EXONhg19.bed', 'r').readlines()

random.shuffle(bedLines)
lineN = 20 if len(sys.argv) < 2 else int(sys.argv[1])

print(''.join(bedLines[:lineN]))

nextflow文件内容如下：

params.bowtie =  "/Users/test/bowtie-1.2.2-macos-x86_64/bowtie"
params.btindex = "/Users/test/Reference/hg19_bt1Index"

params.fq   = "test_1.fq"
params.bed  = "test1.bed"

bed = file(params.bed)
index = file(params.btindex)
fq_ch = Channel.fromPath(params.fq)

process align2BAM {
    input:
    file input_fq from fq_ch

    output:
    file "sorted.bam" into bamfile

    """
    ${params.bowtie} -p 8 -l 5 --chunkmbs 512 -v 3 -k 10 $index $input_fq -S out_1.sam 1> align.log 2>&1
    samtools view -bS out_1.sam | samtools sort - -T /tmp/sortedS -o sorted.bam
    """
}

process plotCovStats {
    input:
    file input_bam from bamfile

    output:
    file "out_plot-quals.png" into result_pic

    """
    samtools stats -t $bed $input_bam > bam_results.stats
    perl plot_bamstats.pl -p out_plot bam_results.stats
    """
}

result.println { it.trim() }

最后使用nextflow的可执行程序 运行该nf文件

./nextflow run tutorial.nf

plot_bamstats.pl 来自以下链接 (实际需要把该perl脚本的绝对路径定义为一个参数， 因为nextflow运行时是在一个新建的work文件夹中运行的)

https://github.com/samtools/samtools/blob/c18f0f59c84b8c343420713547984aec1709729f/misc/plot-bamstats

最终会有一个 结果图片out_plot-quals.png, 如果没有安装gnuplot，则无法正常运行。