ChIP-seq生物信息学科研信息学

ChIP-MACS2大批量运行

2019-04-03  本文已影响16人  caokai001

当需要进行几百个peak calling 时候。对对应的treat,control bam写到macs2_input.txt文件中,就可以批量跑。
类似这种三列形式;
“输出结果文件名”,“treat.bam” ,“control.bam”


image.png

但是文件太多,可能文件名输入错误,所以需要检测一下macs2_input.txt中文件是否都存在,以免出错。

[kcao@login 3.bam]$ cat Macs2_uniq.bam.txt|while read id;
do A=($id);for i in `seq 1 $[ ${#A[*]}-1]`;
do [ -f ${A[i]} ] && echo "${A[i]}--->found" || echo "${A[i]}--->not found";
done;
done|grep "not found"

LNCaP_input_H3k27ac_DHT.uniq.bam--->not found
LNCaP_input_H3k27ac_DHT.uniq.bam--->not found

呀!发现了此文件找不到,仔细找找原因~~

下面就是批量跑macs2
ps:和上面代码没衔接,理解一下就好

cat Macs2_list.txt|while read id;
do
    echo $id
    arr=($id);
    treat=${arr[1]};
    control=${arr[2]};
    name=${arr[0]};
    if [ -z $control ];then 
        echo "macs2 callpeak -t ../${treat}.uniq.bam -g hs -n ./1.narrow/${name}_narrow"|qsub -d ./ -N ${name}.narrow
        echo "macs2 callpeak -t ../${treat}.uniq.bam --broad -g hs --broad-cutoff 0.1 -n ./2.broad/${name}_broad"|qsub -d ./ -N ${name}.broad
    else
        echo "macs2 callpeak -t ../${treat}.uniq.bam -c ../${control}.uniq.bam -g hs -n ./1.narrow/${name}_narrow"|qsub -d ./ -N ${name}.narrow
        echo "macs2 callpeak -t ../${treat}.uniq.bam -c ../${control}.uniq.bam -g hs -n ./2.broad/${name}_broad --broad --broad-cutoff 0.1 "|qsub -d ./ -N ${name}.broad
    fi
done

ok~

上一篇 下一篇

猜你喜欢

热点阅读