ticket_296259/run.sh at main · UPPMAX/ticket_296259 · 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
129
130
131
132
133
134
135
136
137
138
139
140
#!/bin/bash -l
#SBATCH -A staff
#SBATCH -t 1:00:00
#SBATCH -p devcore

module load bioinfo-tools
module load samtools
module load vartrix
module load java/OpenJDK_17+35

# Added by Richel
module load picard

# Added by Richel
echo "Skip the first part"

# Added by Richel
if [[ 1 == 2 ]]; then

    # Commented out by Richel
    # <<'_SKIP_'

    # Original bam
    cp ../8_EXPRESSION.D/output.d/01_out/05.0_rsem/SS2_19_037/H13/SS2_19_037-H13.genome.sorted.bam ./H13/SS2_19_037-H13.genome.sorted-original.bam
    cp ../8_EXPRESSION.D/output.d/01_out/05.0_rsem/SS2_19_037/H13/SS2_19_037-H13.genome.sorted.bam.bai .H13/SS2_19_037-H13.genome.sorted-original.bam
    bamO="SS2_19_037-H13.genome.sorted-original.bam"
    bai="SS2_19_037-H13.genome.sorted-original.bam.bai"

    # Split to chr2
    ## aligment
    samtools view -b ${bamO} chr2 > SS2_19_037-H13_chr2.bam


    ## fasta reference
    fastafile="/castor/project/proj/maria.d/8_EXPRESSION.D/09.0_vartrix/data-d/ref_fasta-d/00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92.fa"
    cp ${fastafile} .
    samtools faidx 00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92.fa chr2 >  ./00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92-chr2.fa
    samtools faidx ./00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92-chr2.fa
    rm 00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92.fa.fai

    # Commented out by Richel
    # _SKIP_
fi

# Added by Richel
echo "Add missing tag"

# Add missing tag
tooldir="/home/mararc/bin/jvarkit/dist/"
LINE="./SS2_19_037-H13_chr2.bam"
outdir1="./"
cell="SS2_19_037-H13"

# Added by Richel
echo "Running jvarkit"
if [ -f ${tooldir}jvarkit.jar ]; then

    # Added by Richel
    echo "Running jvarkit in user's way"

    java -jar ${tooldir}jvarkit.jar samjdk -e 'String c=record.getReadName(); int h=0; int s=21; record.setAttribute("CB",c.substring(h,s));return record;' ${LINE}  > ${outdir1}${cell}-CB_chr2.sam
else

    # Added by Richel
    echo "Check if jvarkit can be run in Richel's way"

    if [ ! -f jvarkit.sif ]; then
        echo "ERROR: 'jvarkit.sif' not found. See https://docs.uppmax.uu.se/software/jvarkit how to create it"
        exit 42
    fi

    # Added by Richel
    echo "Running jvarkit in Richel's way"

    ./jvarkit.sif java -jar /opt/jvarkit/dist/jvarkit.jar samjdk -e 'String c=record.getReadName(); int h=0; int s=21; record.setAttribute(\"CB\",c.substring(h,s));return record;' ${LINE}  > ${outdir1}${cell}-CB_chr2.sam

    # Added by Richel
    echo "Done running jvarkit in Richel's way"

fi

# Added by Richel
echo "Convert sam to bam"


# convert back from sam to bam -- I met with Richel @UPPMAX-SUPPORT and he pointed out the file jvarkit spits out is a text, he's right. I actually get a sam format and not bam
# Therefore, I will now convert back
samtools view -bS ./SS2_19_037-H13-CB_chr2.sam > ./SS2_19_037-H13_chr2_bam.bam

# Added by Richel
echo "Indexing"

## and index
samtools index SS2_19_037-H13_chr2_bam.bam

# Added by Richel
echo "Done indexing"

################# Added by Richel, start
echo "----------------------------------------------"
echo "Fix problem here"
# The file, SS2_19_037-H13-CB_chr2.sam, is already ready:
head SS2_19_037-H13-CB_chr2.sam


echo "----------------------------------------------"
################# Added by Richel, end


# Added by Richel
echo "Setting vcffile"

# Run - it works!! The only thing left to do is add the tag to the bam header sth like @CB:ZXXXXX
# reference vcf file
vcffile="./ALK_chr.vcf"

# Added by Richel
echo "Setting fastafile"

# fasta file for the reference genome used
fastafile="./00.0_chrom_seq_removed_GRCh38.primary_assembly.genome-nochrY_ERCC92-chr2.fa"

# Added by Richel
echo "Setting LINE"

LINE="./"

# Added by Richel
echo "Setting cell"

cell="SS2_19_037-H13"

# Added by Richel
echo "Running VarTrix"

vartrix --vcf ${vcffile} --bam ${LINE}${cell}_chr2_bam.bam --fasta ${fastafile} --cell-barcodes ${LINE}barcodes --out-matrix ./out_mat --out-variants ./out_var &> log


samtools view -S -b sample.sam > sample.bam