MURVseq/enzymeseq_pipe_13_1.sh at main · chorong5155/MURVseq · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#!/bin/bash
topdir=$(pwd)
account="p31754"
queue="long"

OUT_DIR="$topdir/natural_TSS"
natural="$topdir/un_natural_output"
GENOME="/projects/b1154/referenceGenomes/hg38+VacV/hg38+VacV"
BED_FILE="/projects/b1154/referenceGenomes/hg38+VacV/VacV_150utr_flanks.bed"
GENOME_FA="/projects/b1154/referenceGenomes/hg38+VacV/hg38+VacV.fa"

module load bedtools bowtie2 samtools
mkdir -p "$OUT_DIR" "$topdir/debug"

for sample_dir in "$natural"/*/; do
    sample=$(basename "$sample_dir")
    sample_out="$OUT_DIR/$sample"
    mkdir -p "$sample_out"

    sbatch --export=ALL,SAMPLE="$sample",SAMPLE_DIR="$sample_dir",OUT_DIR="$sample_out",GENOME="$GENOME",BED_FILE="$BED_FILE",GENOME_FA="$GENOME_FA" <<EOF
#!/bin/bash -l
#SBATCH -A $account
#SBATCH -p $queue
#SBATCH -o $topdir/debug/TSS_${sample}.out
#SBATCH -e $topdir/debug/TSS_${sample}.err
#SBATCH -t 24:00:00
#SBATCH -c 2
#SBATCH --mem-per-cpu=40G
#SBATCH -J TSS_${sample}

module load bedtools samtools minimap2

EMPTY_REPORT="\$OUT_DIR/empty_trimmed_counts.tsv"

for FASTA in "\$SAMPLE_DIR"/naturalUTRs/*_naturalUTRs.fasta; do
    base=\$(basename "\$FASTA" .fasta)
    PREFIX="\$OUT_DIR/\$base"


    bowtie2 --local -L 14 -f -x "\$GENOME" -U "\${FASTA}" -S "\${PREFIX}.sam"

    samtools view -bS "\${PREFIX}.sam" | samtools sort -o "\${PREFIX}.sorted.bam"
    samtools index "\${PREFIX}.sorted.bam"
    samtools view "\${PREFIX}.sorted.bam" |
    awk '
        # Skip unmapped reads
        (\$2 % 8) >= 4 { next }

        {
            chr=\$3; pos=\$4; read=\$1; flag=\$2; cigar=\$6

            # compute reference-consuming length
            refLen=0
            cig=cigar
            while (match(cig, /[0-9]+[MIDNSHP=X]/)) {
                c=substr(cig, RSTART, RLENGTH)
                gsub(/[A-Z=]/, "", c)
                op=substr(cig, RSTART+length(c), 1)
                n=c+0

                if (op=="M" || op=="D" || op=="N" || op=="=" || op=="X")
                    refLen += n

                cig=substr(cig, RSTART+RLENGTH)
            }

            # strand logic (poxvirus)
            if (flag % 32 >= 16) {
                strand = "-"
                tss = pos + refLen - 1
            } else {
                strand = "+"
                tss = pos
            }

            start = tss - 1
            if (start < 0) start = 0

            print chr "\t" start "\t" tss "\t" read "\t.\t" strand
        }
    ' > "\${PREFIX}_TSS.bed"

awk -v OFS="\t" '{ key=\$1"\t"\$2"\t"\$3; counts[key]++ }
END {
    for (k in counts) print k, counts[k]
}' "\${PREFIX}_TSS.bed" \
| sort -k1,1 -k2,2n > "\${PREFIX}_TSS_unique.bed"


    bedtools intersect -a "\${PREFIX}_TSS.bed" -b "\$BED_FILE" -wa -wb \
        > "\${PREFIX}_TSS_annotated.bed"

    #cut -f1-3 "\${PREFIX}_TSS.bed" | sort -u > "\${PREFIX}_TSS_unique.bed"

awk -v OFS="\t" '
{
    key=\$1"\t"\$3"\t"\$6
    counts[key]++
}
END {
    for (k in counts) {
        print counts[k], k
    }
}
' "\${PREFIX}_TSS.bed" \
| sort -nrk1,1 | head -3 \
| awk -v OFS="\t" '{print \$2, \$3, \$4, \$1}' \
> "\${PREFIX}_dominant_TSS.txt"


    bedtools getfasta \
        -fi "\${GENOME_FA}" \
        -bed "\${PREFIX}_TSS.bed" \
        -s \
        -name \
        -fo "\${PREFIX}_TSS_seq.fa"

    awk -v OFS="\t" -v stage="early" '{
        TSS_ID="TSS_"NR
        print \$1, \$2, \$3, TSS_ID, ".", \$6, \$7, \$8, \$9, \$10, stage
    }' "\${PREFIX}_TSS_annotated.bed" > "\${PREFIX}_TSS_final.tsv"


done
EOF
done