Skip to content

Commit 31e2084

Browse files
committed
added synthetic fastqs for fusion detection
1 parent 7dcb3dd commit 31e2084

3 files changed

Lines changed: 102 additions & 0 deletions

File tree

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
#!/usr/bin/env bash
2+
set -euo pipefail
3+
4+
FASTA="minigenome.fa"
5+
OUTDIR="synthetic_reads"
6+
7+
mkdir -p "$OUTDIR"/{genes,fusions}
8+
9+
echo "====================================="
10+
echo "Step 1: index FASTA"
11+
echo "====================================="
12+
samtools faidx "$FASTA"
13+
14+
echo "====================================="
15+
echo "Step 2: define gene list (NO coordinates)"
16+
echo "====================================="
17+
18+
GENES=(
19+
AKAP9 ALK ATF1 BRAF BRD4 CD74 EML4 ETV1 ETV6
20+
EWSR1 FLI1 HOOK3 NTRK3 NUTM1 RET ROS1 TMPRSS2
21+
)
22+
23+
echo "====================================="
24+
echo "Step 3: build normal gene transcripts"
25+
echo "====================================="
26+
27+
for gene in "${GENES[@]}"; do
28+
seq=$(samtools faidx "$FASTA" "$gene" | tail -n +2)
29+
30+
echo ">${gene}_normal" > "$OUTDIR/genes/${gene}.fa"
31+
echo "$seq" >> "$OUTDIR/genes/${gene}.fa"
32+
done
33+
34+
cat "$OUTDIR"/genes/*.fa > "$OUTDIR/all_genes.fa"
35+
36+
echo "====================================="
37+
echo "Step 4: build fusion transcripts"
38+
echo "====================================="
39+
40+
FUSIONS=(
41+
"ALK EML4"
42+
"ROS1 CD74"
43+
"RET NTRK3"
44+
"BRAF ETV6"
45+
"EWSR1 FLI1"
46+
"TMPRSS2 ALK"
47+
)
48+
49+
i=0
50+
for pair in "${FUSIONS[@]}"; do
51+
g1=$(echo $pair | awk '{print $1}')
52+
g2=$(echo $pair | awk '{print $2}')
53+
54+
seq1=$(samtools faidx "$FASTA" "$g1" | tail -n +2)
55+
seq2=$(samtools faidx "$FASTA" "$g2" | tail -n +2)
56+
57+
mid1=$(( ${#seq1} / 2 ))
58+
mid2=$(( ${#seq2} / 3 ))
59+
60+
fusion_seq="${seq1:0:$mid1}${seq2:$mid2}"
61+
62+
echo ">fusion_${g1}_${g2}_${i}" > "$OUTDIR/fusions/fusion_${i}.fa"
63+
echo "$fusion_seq" >> "$OUTDIR/fusions/fusion_${i}.fa"
64+
65+
i=$((i+1))
66+
done
67+
68+
cat "$OUTDIR"/fusions/*.fa > "$OUTDIR/all_fusions.fa"
69+
70+
echo "====================================="
71+
echo "Step 5: combine transcriptome"
72+
echo "====================================="
73+
74+
cat "$OUTDIR/all_genes.fa" "$OUTDIR/all_fusions.fa" > "$OUTDIR/transcriptome.fa"
75+
76+
echo "====================================="
77+
echo "Step 6: simulate reads"
78+
echo "====================================="
79+
80+
art_illumina \
81+
-ss HS25 \
82+
-i "$OUTDIR/all_genes.fa" \
83+
-p -l 150 -f 30 -m 200 -s 10 -rs 42 \
84+
-o "$OUTDIR/normal_"
85+
86+
art_illumina \
87+
-ss HS25 \
88+
-i "$OUTDIR/all_fusions.fa" \
89+
-p -l 150 -f 5 -m 200 -s 10 -rs 42 \
90+
-o "$OUTDIR/fusion_"
91+
92+
echo "====================================="
93+
echo "Step 7: merge reads"
94+
echo "====================================="
95+
96+
cat "$OUTDIR/normal_1.fq" "$OUTDIR/fusion_1.fq" > "$OUTDIR/final_1.fq"
97+
cat "$OUTDIR/normal_2.fq" "$OUTDIR/fusion_2.fq" > "$OUTDIR/final_2.fq"
98+
99+
gzip -f "$OUTDIR/final_1.fq"
100+
gzip -f "$OUTDIR/final_2.fq"
101+
102+
echo "DONE"

fastqs/test_sample1_R1.fastq.gz

44 MB
Binary file not shown.

fastqs/test_sample1_R2.fastq.gz

45.3 MB
Binary file not shown.

0 commit comments

Comments
 (0)