2020-07-16 靶向捕获测序数据分析记录3

昨天一直在调试代码，之前学习的时候只有一个样本，运行时间比较短，没有遇到现在的问题，现在实际分析数据是有600多个样本的，平均一个样本25分钟的样子，仅转换版本号可能预计花费时间约12天。所以必须提交到后台运行。

我之前的基础学习并没有深入了解nohup指令，直接把它写入循环了，而且没有写判断语句，因此出了不少的错。

这里总结一下，首先，完整的代码应该是nohup command &,&不能少。

所以应该把循环写入脚本然后用nohup提交到后台运行

版本转换脚本crossmap.sh
用vi编辑，然后chmod 777 crossmap.sh修改属性。

##写成判断语句
cat config | while read id
do 
    bam=~/CHD_pooling_seq/${id}.dedup.bam
    if [ ! -f ~/project/0.bwa/ok.${id}_marked.status ]
    then
        echo "start CrossMap for ${id}" `date`
        python /root/miniconda3/envs/py3/bin/CrossMap.py \
        bam ~/biosoft/liftover/hg19ToHg38.over.chain.gz \
        ${bam} \
        ~/project/0.bwa/${id}.hg38 \
        1>~/project/0.bwa/${id}_log.mark 2>&1
        
        if [ $? -eq 0 ]
        then
            touch ~/project/0.bwa/ok.${id}_marked.status
        fi
        echo "end CrossMap for ${id}" `date`
     fi 
done

然后再nohup bash crossmap.sh &将任务提交到后台。

(base) root@1100150:~/project# ll
total 1252
drwxr-xr-x  8 root root   193 Jul 15 10:07 ./
drwx------ 15 root root    21 Jul 13 11:41 ../
drwxr-xr-x  2 root root   789 Jul 16 13:21 0.bwa/
drwxr-xr-x  4 root root     4 Jul 13 06:22 1.qc/
drwxr-xr-x  3 root root  1383 Jul 15 10:07 2.gatk/
drwxr-xr-x  6 root root     6 Jul 10 00:43 3.annotation/
-rwxrwxrwx  1 root root  1170 Jul 15 08:22 bqsr.sh*
-rw-r--r--  1 root root  2937 Jul 13 05:44 config
-rwxrwxrwx  1 root root   538 Jul 15 08:11 crossmap.sh*

image.png

运行2天后的结果如下：

(base) root@1100150:~# cd ~/project/0.bwa
(base) root@1100150:~/project/0.bwa# l |grep 'ok'
C163_log.mark              ok.C100_marked.status
C164.hg38.bam              ok.C102_marked.status
C164.hg38.sorted.bam           ok.C103_marked.status
C164.hg38.sorted.bam.bai       ok.C104_marked.status
C164_log.mark              ok.C106_marked.status
C165.hg38.bam              ok.C107_marked.status
C165.hg38.sorted.bam           ok.C108_marked.status
C165.hg38.sorted.bam.bai       ok.C109_marked.status
C165_log.mark              ok.C10_marked.status
C166.hg38.bam              ok.C110_marked.status
C166.hg38.sorted.bam           ok.C111_marked.status
C166.hg38.sorted.bam.bai       ok.C112_marked.status
C166_log.mark              ok.C113_marked.status
C167.hg38.bam              ok.C114_marked.status
C167.hg38.sorted.bam           ok.C115_marked.status
C167.hg38.sorted.bam.bai       ok.C116_marked.status
C167_log.mark              ok.C117_marked.status
C168.hg38.bam              ok.C118_marked.status
C168.hg38.sorted.bam           ok.C119_marked.status
C168.hg38.sorted.bam.bai       ok.C120_marked.status
C168_log.mark              ok.C122_marked.status
C169.hg38.bam              ok.C124_marked.status
C169.hg38.sorted.bam           ok.C126_marked.status
C169.hg38.sorted.bam.bai       ok.C127_marked.status
C169_log.mark              ok.C128_marked.status
C16_log.mark               ok.C129_marked.status
C17.hg38.bam               ok.C12_marked.status
C17.hg38.sorted.bam        ok.C130_marked.status
C17.hg38.sorted.bam.bai        ok.C131_marked.status
C170.hg38.bam              ok.C132_marked.status
C170.hg38.sorted.bam           ok.C135_marked.status
C170.hg38.sorted.bam.bai       ok.C136_marked.status
C170_log.mark              ok.C13_marked.status
C171.hg38.bam              ok.C140_marked.status
C171.hg38.sorted.bam           ok.C141_marked.status
C171.hg38.sorted.bam.bai       ok.C142_marked.status
C171_log.mark              ok.C143_marked.status
C173.hg38.bam              ok.C144_marked.status
C173.hg38.sorted.bam.tmp.0000.bam  ok.C145_marked.status
C173_log.mark              ok.C146_marked.status
C174.hg38.bam              ok.C147_marked.status
C174.hg38.sorted.bam           ok.C148_marked.status
C174.hg38.sorted.bam.bai       ok.C149_marked.status
C174_log.mark              ok.C14_marked.status
C175.hg38.bam              ok.C150_marked.status
C175.hg38.sorted.bam           ok.C152_marked.status
C175.hg38.sorted.bam.bai       ok.C153_marked.status
C175_log.mark              ok.C155_marked.status
C176.hg38.bam              ok.C156_marked.status
C176.hg38.sorted.bam           ok.C157_marked.status
C176.hg38.sorted.bam.bai       ok.C158_marked.status
C176_log.mark              ok.C159_marked.status
C177.hg38.bam              ok.C15_marked.status
C177.hg38.sorted.bam           ok.C160_marked.status
C177.hg38.sorted.bam.bai       ok.C161_marked.status
C177_log.mark              ok.C162_marked.status
C17_log.mark               ok.C163_marked.status
C180.hg38.bam              ok.C164_marked.status
C180.hg38.sorted.bam           ok.C165_marked.status
C180.hg38.sorted.bam.bai       ok.C166_marked.status
C180_log.mark              ok.C167_marked.status
C181.hg38.bam              ok.C168_marked.status
C181.hg38.sorted.bam           ok.C169_marked.status
C181.hg38.sorted.bam.bai       ok.C16_marked.status
C181_log.mark              ok.C170_marked.status
C182.hg38.bam              ok.C171_marked.status
C182.hg38.sorted.bam           ok.C173_marked.status
C182.hg38.sorted.bam.bai       ok.C174_marked.status
C182_log.mark              ok.C175_marked.status
C183.hg38.bam              ok.C176_marked.status
C183.hg38.sorted.bam           ok.C177_marked.status
C183.hg38.sorted.bam.bai       ok.C17_marked.status
C183_log.mark              ok.C180_marked.status
C184.hg38.bam              ok.C181_marked.status
C184.hg38.sorted.bam           ok.C182_marked.status
C184.hg38.sorted.bam.bai       ok.C183_marked.status
C184_log.mark              ok.C184_marked.status
C185.hg38.bam              ok.C185_marked.status
C185.hg38.sorted.bam           ok.C186_marked.status
C185.hg38.sorted.bam.bai       ok.C187_marked.status
C185_log.mark              ok.C188_marked.status
C186.hg38.bam              ok.C190_marked.status
C186.hg38.sorted.bam           ok.C191_marked.status
C186.hg38.sorted.bam.bai       ok.C192_marked.status
C186_log.mark              ok.C193_marked.status
C187.hg38.bam              ok.C194_marked.status
C187.hg38.sorted.bam           ok.C195_marked.status
C187.hg38.sorted.bam.bai       ok.C196_marked.status
C187_log.mark              ok.C197_marked.status
C188.hg38.bam              ok.C198_marked.status
C188.hg38.sorted.bam           ok.C199_marked.status
C188.hg38.sorted.bam.bai       ok.C19_marked.status
C188_log.mark              ok.C1_marked.status
C19.hg38.bam               ok.C200_marked.status
C19.hg38.sorted.bam        ok.C201_marked.status
C19.hg38.sorted.bam.bai        ok.C202_marked.status
C190.hg38.bam              ok.C203_marked.status
C190.hg38.sorted.bam           ok.C204_marked.status
C190.hg38.sorted.bam.bai       ok.C205_marked.status
C190_log.mark              ok.C206_marked.status
C191.hg38.bam              ok.C207_marked.status
C191.hg38.sorted.bam           ok.C208_marked.status
C191.hg38.sorted.bam.bai       ok.C209_marked.status
C191_log.mark              ok.C210_marked.status
C192.hg38.bam              ok.C211_marked.status
C192.hg38.sorted.bam           ok.C212_marked.status
C192.hg38.sorted.bam.bai       ok.C213_marked.status
C192_log.mark              ok.C214_marked.status
C193.hg38.bam              ok.C215_marked.status
C193.hg38.sorted.bam           ok.C216_marked.status
C193.hg38.sorted.bam.bai       ok.C217_marked.status
C193_log.mark              ok.C218_marked.status
C194.hg38.bam              ok.C219_marked.status
C194.hg38.sorted.bam           ok.C21_marked.status
C194.hg38.sorted.bam.bai       ok.C220_marked.status
C194_log.mark              ok.C221_marked.status
C195.hg38.bam              ok.C222_marked.status
C195.hg38.sorted.bam           ok.C223_marked.status
C195.hg38.sorted.bam.bai       ok.C224_marked.status
C195_log.mark              ok.C225_marked.status
C196.hg38.bam              ok.C226_marked.status
C196.hg38.sorted.bam           ok.C227_marked.status
C196.hg38.sorted.bam.bai       ok.C228_marked.status
C196_log.mark              ok.C229_marked.status
C197.hg38.bam              ok.C22_marked.status
C197.hg38.sorted.bam           ok.C233_marked.status
C197.hg38.sorted.bam.bai       ok.C234_marked.status
C197_log.mark              ok.C235_marked.status
C198.hg38.bam              ok.C236_marked.status
C198.hg38.sorted.bam           ok.C237_marked.status
C198.hg38.sorted.bam.bai       ok.C238_marked.status
C198_log.mark              ok.C239_marked.status
C199.hg38.bam              ok.C23_marked.status
C199.hg38.sorted.bam           ok.C240_marked.status
C199.hg38.sorted.bam.bai       ok.C24_marked.status
C199_log.mark              ok.C25_marked.status
C19_log.mark               ok.C26_marked.status
C1_log.mark            ok.C27_marked.status
C2.hg38.bam            ok.C28_marked.status
C2.hg38.sorted.bam         ok.C29_marked.status
C2.hg38.sorted.bam.bai         ok.C2_marked.status
C200.hg38.bam              ok.C30_marked.status
C200.hg38.sorted.bam           ok.C31_marked.status
C200.hg38.sorted.bam.bai       ok.C32_marked.status
C200_log.mark              ok.C33_marked.status
C201.hg38.bam              ok.C34_marked.status
C201.hg38.sorted.bam           ok.C35_marked.status
C201.hg38.sorted.bam.bai       ok.C36_marked.status
C201_log.mark              ok.C37_marked.status
C202.hg38.bam              ok.C38_marked.status
C202.hg38.sorted.bam           ok.C39_marked.status
C202.hg38.sorted.bam.bai       ok.C3_marked.status
C202_log.mark              ok.C40_marked.status
C203.hg38.bam              ok.C41_marked.status
C203.hg38.sorted.bam           ok.C42_marked.status
C203.hg38.sorted.bam.bai       ok.C43_marked.status
C203_log.mark              ok.C44_marked.status
C204.hg38.bam              ok.C4_marked.status

同样的方法进行碱基质量值矫正

(wes) root@1100150:~/project# cat bqsr.sh
##已经去重了，下一步进行质量值矫正
GATK=~/biosoft/gatk-4.1.7.0/gatk
#references
ref=~/reference/genome/Homo_sapiens_assembly38.fasta
gatk_ref=~/reference/genome/Homo_sapiens_assembly38.fasta
gatk_bundle=~/annotation/variation/GATK
dbsnp=$gatk_bundle/dbsnp_146.hg38.vcf.gz
indel=$gatk_bundle/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
G1000=$gatk_bundle/1000G_phase1.snps.high_confidence.hg38.vcf.gz
hapmap=$gatk_bundle/hapmap_3.3.hg38.vcf.gz
omini=$gatk_bundle/1000G_omni2.5.hg38.vcf.gz

cat config  | while read id
do
    if [ ! -f ./2.gatk/${id}_bqsr.bam ]
    then
        echo "start BQSR for ${id}" `date`
        $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./"  BaseRecalibrator \
        -R $ref  \
        -I ./0.bwa/${id}.hg38.sorted.bam  \
        --known-sites ${dbsnp} \
        --known-sites ${indel} \
        --known-sites ${G1000} \
        -O ./2.gatk/${id}_recal.table \
        1>./2.gatk/${id}_log.recal 2>&1

        $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./"  ApplyBQSR \
        -R $ref  \
        -I ./0.bwa/${id}.hg38.sorted.bam  \
        -bqsr ./2.gatk/${id}_recal.table \
        -O ./2.gatk/${id}_bqsr.bam \
        1>./2.gatk/${id}_log.ApplyBQSR  2>&1

        echo "end BQSR for ${id}" `date`
    fi
done

查看矫正后的bam的头文件：

(wes) root@1100150:~/project# samtools view -h 2.gatk/C141_bqsr.bam |less -S

image.png