1234567891011121314151617181920212223242526272829303132333435363738 |
- r1_trim=$1
- r2_trim=$2
- sample=$3
- outdir=$4
- curPath=$(dirname $(readlink -f "$0"))
- subPath=$(dirname $(readlink -f $curPath))
- no_IPC_target=${subPath}/DB/no_IPC_target.bed
- IPC=${subPath}/DB/IPC.bed
- ref=/cgdata/IVDRD/18gene_DB/qaseq_fusion/hg38.fa
- mkdir -p ${outdir}/2.mapping;
- bwa mem -t 16 ${ref} ${r1_trim} ${r2_trim} -R "@RG\tID:1\tLB:lib1\tPL:ILLUMINA\tSM:${sample}\tPU:unit1" > ${outdir}/2.mapping/${sample}.sam;
- samtools view ${outdir}/2.mapping/${sample}.sam --threads 8 -q 30 -bS > ${outdir}/2.mapping/${sample}.pe.bam;
- samtools view ${outdir}/2.mapping/${sample}.pe.bam > ${outdir}/2.mapping/${sample}.q30.sam;
- sambamba sort ${outdir}/2.mapping/${sample}.pe.bam --tmpdir=${outdir}/tmp -o ${outdir}/2.mapping/${sample}.sort.bam -t 8;
- samtools view -L $no_IPC_target ${outdir}/2.mapping/${sample}.sort.bam > ${outdir}/2.mapping/${sample}.sort.target.sam;
- samtools view -L $IPC ${outdir}/2.mapping/${sample}.sort.bam > ${outdir}/2.mapping/${sample}.sort.IPC.sam;
- awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.target.sam | sort | uniq > ${outdir}/2.mapping/${sample}.target.reads.txt;
- awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.IPC.sam | sort | uniq > ${outdir}/2.mapping/${sample}.IPC.reads.txt;
- awk 'NR==FNR{a[$1]="OK";next}{print $0"\t"a[$1]}' ${outdir}/2.mapping/${sample}.IPC.reads.txt ${outdir}/2.mapping/${sample}.q30.sam | awk '{if($NF!="OK"){print $0}}' > ${outdir}/2.mapping/${sample}.q30.tmp.sam;
- awk 'NR==FNR{a[$1]="OK";next}{print $0"\t"a[$1]}' ${outdir}/2.mapping/${sample}.target.reads.txt ${outdir}/2.mapping/${sample}.q30.tmp.sam | awk '{if($NF!="OK"){print $0}}' > ${outdir}/2.mapping/${sample}.sort.nonspecific.sam;
- awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.nonspecific.sam | sort | uniq > ${outdir}/2.mapping/${sample}.nonspecific.reads.txt
|