Skip to content

Commit 1d96136

Browse files
khersameesh24khersameesh24
committed
repo cleanup
1 parent 1ce4f8f commit 1d96136

4 files changed

Lines changed: 234 additions & 51 deletions

File tree

‎.gitignore‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
ref/Homo_sapiens.GRCh38.110.gtf
2+
ref/Homo_sapiens.GRCh38.dna.primary_assembly.fa
3+
ref/Homo_sapiens.GRCh38.dna.primary_assembly.fa.fai

‎fastqs/generate_synthetic_fastqs.sh‎

Lines changed: 82 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -6,36 +6,21 @@ OUTDIR="synthetic_reads"
66

77
mkdir -p "$OUTDIR"/{genes,fusions}
88

9-
echo "====================================="
10-
echo "Step 1: index FASTA"
11-
echo "====================================="
9+
echo "== Index FASTA =="
1210
samtools faidx "$FASTA"
1311

14-
echo "====================================="
15-
echo "Step 2: define gene list (NO coordinates)"
16-
echo "====================================="
17-
1812
GENES=(
1913
AKAP9 ALK ATF1 BRAF BRD4 CD74 EML4 ETV1 ETV6
2014
EWSR1 FLI1 HOOK3 NTRK3 NUTM1 RET ROS1 TMPRSS2
2115
)
2216

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"
17+
extract_gene() {
18+
samtools faidx "$FASTA" "$1" 2>/dev/null \
19+
| tail -n +2 \
20+
| tr -d '\n'
21+
}
3522

36-
echo "====================================="
37-
echo "Step 4: build fusion transcripts"
38-
echo "====================================="
23+
echo "== Build breakpoint fusion fragments =="
3924

4025
FUSIONS=(
4126
"ALK EML4"
@@ -46,57 +31,103 @@ FUSIONS=(
4631
"TMPRSS2 ALK"
4732
)
4833

34+
rm -f "$OUTDIR"/fusions/*.fa
35+
4936
i=0
37+
5038
for pair in "${FUSIONS[@]}"; do
51-
g1=$(echo $pair | awk '{print $1}')
52-
g2=$(echo $pair | awk '{print $2}')
5339

54-
seq1=$(samtools faidx "$FASTA" "$g1" | tail -n +2)
55-
seq2=$(samtools faidx "$FASTA" "$g2" | tail -n +2)
40+
g1=$(echo "$pair" | awk '{print $1}')
41+
g2=$(echo "$pair" | awk '{print $2}')
42+
43+
seq1=$(extract_gene "$g1")
44+
seq2=$(extract_gene "$g2")
45+
46+
[[ -z "$seq1" || -z "$seq2" ]] && continue
47+
48+
bp1=$(( ${#seq1} / 2 ))
49+
bp2=$(( ${#seq2} / 3 ))
5650

57-
mid1=$(( ${#seq1} / 2 ))
58-
mid2=$(( ${#seq2} / 3 ))
51+
# 200 bp on each side of breakpoint
52+
left_start=$(( bp1 - 200 ))
53+
(( left_start < 0 )) && left_start=0
5954

60-
fusion_seq="${seq1:0:$mid1}${seq2:$mid2}"
55+
left_seq="${seq1:$left_start:200}"
56+
right_seq="${seq2:$bp2:200}"
6157

62-
echo ">fusion_${g1}_${g2}_${i}" > "$OUTDIR/fusions/fusion_${i}.fa"
63-
echo "$fusion_seq" >> "$OUTDIR/fusions/fusion_${i}.fa"
58+
fusion_fragment="${left_seq}${right_seq}"
6459

65-
i=$((i+1))
60+
{
61+
echo ">fusion_${g1}_${g2}_${i}"
62+
echo "$fusion_fragment"
63+
} > "$OUTDIR/fusions/fusion_${i}.fa"
64+
65+
((i+=1))
6666
done
6767

6868
cat "$OUTDIR"/fusions/*.fa > "$OUTDIR/all_fusions.fa"
6969

70-
echo "====================================="
71-
echo "Step 5: combine transcriptome"
72-
echo "====================================="
73-
74-
cat "$OUTDIR/all_genes.fa" "$OUTDIR/all_fusions.fa" > "$OUTDIR/transcriptome.fa"
70+
echo "== Simulate normal background reads =="
7571

76-
echo "====================================="
77-
echo "Step 6: simulate reads"
78-
echo "====================================="
72+
# Keep background modest.
73+
# With a 4.4 Mb minigenome this produces ~40-60 MB gzipped FASTQs.
7974

8075
art_illumina \
8176
-ss HS25 \
82-
-i "$OUTDIR/all_genes.fa" \
83-
-p -l 150 -f 30 -m 200 -s 10 -rs 42 \
77+
-i "$FASTA" \
78+
-p \
79+
-l 150 \
80+
-f 50 \
81+
-m 250 \
82+
-s 25 \
83+
-rs 42 \
8484
-o "$OUTDIR/normal_"
8585

86+
echo "== Simulate breakpoint-spanning fusion reads =="
87+
88+
# High coverage on tiny breakpoint fragments
89+
# Generates many split/chimeric reads without huge files
90+
8691
art_illumina \
8792
-ss HS25 \
8893
-i "$OUTDIR/all_fusions.fa" \
89-
-p -l 150 -f 5 -m 200 -s 10 -rs 42 \
94+
-p \
95+
-l 150 \
96+
-f 600 \
97+
-m 250 \
98+
-s 25 \
99+
-rs 43 \
90100
-o "$OUTDIR/fusion_"
91101

92-
echo "====================================="
93-
echo "Step 7: merge reads"
94-
echo "====================================="
102+
echo "== Merge =="
103+
104+
cat \
105+
"$OUTDIR/normal_1.fq" \
106+
"$OUTDIR/fusion_1.fq" \
107+
> "$OUTDIR/test_sample_R1.fastq"
108+
109+
cat \
110+
"$OUTDIR/normal_2.fq" \
111+
"$OUTDIR/fusion_2.fq" \
112+
> "$OUTDIR/test_sample_R2.fastq"
113+
114+
echo "== Compress =="
115+
116+
gzip -f "$OUTDIR/test_sample_R1.fastq"
117+
gzip -f "$OUTDIR/test_sample_R2.fastq"
118+
119+
echo
120+
echo "Final sizes:"
121+
du -h "$OUTDIR/test_sample_R1.fastq.gz"
122+
du -h "$OUTDIR/test_sample_R2.fastq.gz"
95123

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"
124+
echo
125+
echo "Read counts:"
126+
echo -n "R1: "
127+
zcat "$OUTDIR/test_sample_R1.fastq.gz" | awk 'END{print NR/4}'
98128

99-
gzip -f "$OUTDIR/final_1.fq"
100-
gzip -f "$OUTDIR/final_2.fq"
129+
echo -n "R2: "
130+
zcat "$OUTDIR/test_sample_R2.fastq.gz" | awk 'END{print NR/4}'
101131

102-
echo "DONE"
132+
echo
133+
echo "Done"

‎ref/generate_minigenome.sh‎

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
#!/usr/bin/env bash
2+
set -euo pipefail
3+
4+
### INPUTS
5+
FASTA="Homo_sapiens.GRCh38.dna.primary_assembly.fa"
6+
GTF="Homo_sapiens.GRCh38.110.gtf"
7+
8+
OUTDIR="minigenome_arriba"
9+
THREADS=8
10+
PAD=100000 # window around genes (controls final size)
11+
12+
mkdir -p "$OUTDIR"
13+
14+
echo "=============================="
15+
echo "Step 1: index FASTA"
16+
echo "=============================="
17+
samtools faidx "$FASTA"
18+
19+
echo "=============================="
20+
echo "Step 2: define genes"
21+
echo "=============================="
22+
23+
GENES=(
24+
AKAP9 ALK ATF1 BRAF BRD4 CD74 EML4 ETV1 ETV6
25+
EWSR1 FLI1 HOOK3 NTRK3 NUTM1 RET ROS1 TMPRSS2
26+
)
27+
28+
echo "=============================="
29+
echo "Step 3: build gene BED (FIXED)"
30+
echo "=============================="
31+
32+
printf "%s\n" "${GENES[@]}" > "$OUTDIR/genes.tmp"
33+
34+
awk -F'\t' -v OFS='\t' -v pad="$PAD" '
35+
BEGIN {
36+
while ((getline < "'$OUTDIR'/genes.tmp") > 0) g[$1]=1
37+
}
38+
39+
$3=="gene" {
40+
if (match($9, /gene_name "([^"]+)"/, m)) {
41+
if (m[1] in g) {
42+
start=$4-pad; if (start<1) start=1;
43+
print $1, start, $5+pad, m[1]
44+
}
45+
}
46+
}
47+
' "$GTF" > "$OUTDIR/genes.bed"
48+
49+
rm -f "$OUTDIR/genes.tmp"
50+
51+
echo "=============================="
52+
echo "Step 4: build compact FASTA (FIXED)"
53+
echo "=============================="
54+
55+
> "$OUTDIR/minigenome.fa"
56+
57+
while read -r chr start end gene; do
58+
59+
region="${chr}:${start}-${end}"
60+
echo "Extracting $gene -> $region"
61+
62+
seq=$(samtools faidx "$FASTA" "$region" 2>/dev/null | tail -n +2 || true)
63+
64+
if [[ -z "$seq" ]]; then
65+
echo "WARNING: skipping empty region $region"
66+
continue
67+
fi
68+
69+
{
70+
echo ">$gene"
71+
printf "%s\n" "$seq" | fold -w 60
72+
} >> "$OUTDIR/minigenome.fa"
73+
74+
done < "$OUTDIR/genes.bed"
75+
76+
# HARD VALIDATION (important)
77+
if [[ ! -s "$OUTDIR/minigenome.fa" ]]; then
78+
echo "ERROR: FASTA is empty — check BED extraction"
79+
exit 1
80+
fi
81+
82+
samtools faidx "$OUTDIR/minigenome.fa"
83+
84+
echo "=============================="
85+
echo "Fixing GTF (gene-name based, STAR-safe)"
86+
echo "=============================="
87+
88+
awk -F'\t' -v OFS='\t' '
89+
BEGIN {
90+
while ((getline < "'$OUTDIR'/genes.tmp") > 0) g[$1]=1
91+
}
92+
93+
$3=="exon" {
94+
if (match($9, /gene_name "([^"]+)"/, m)) {
95+
if (m[1] in g) {
96+
# rewrite chromosome to match FASTA contig (gene name)
97+
$1 = m[1]
98+
print
99+
}
100+
}
101+
}
102+
' "$GTF" > "$OUTDIR/minigenome.gtf"
103+
104+
echo "=============================="
105+
echo "Step 6: sanity check (CRITICAL)"
106+
echo "=============================="
107+
108+
echo "FASTA size:"
109+
du -h "$OUTDIR/minigenome.fa"
110+
111+
echo "Contigs:"
112+
cut -f1 "$OUTDIR/minigenome.fa.fai"
113+
114+
echo "GTF features:"
115+
cut -f3 "$OUTDIR/minigenome.gtf" | sort | uniq -c
116+
117+
echo "=============================="
118+
echo "Step 7: STAR index (GitHub-safe settings)"
119+
echo "=============================="
120+
121+
STAR \
122+
--runThreadN "$THREADS" \
123+
--runMode genomeGenerate \
124+
--genomeDir "$OUTDIR/star_index" \
125+
--genomeFastaFiles "$OUTDIR/minigenome.fa" \
126+
--sjdbGTFfile "$OUTDIR/minigenome.gtf" \
127+
--sjdbOverhang 149 \
128+
--genomeSAindexNbases 10 \
129+
--limitGenomeGenerateRAM 20000000000
130+
131+
echo "=============================="
132+
echo "DONE"
133+
echo "=============================="

‎ref/genes.bed‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
ALK 2 29421449 29471450
2+
EML4 2 42503531 42553532
3+
BRAF 7 140462383 140512384
4+
CD74 5 149759242 149809243
5+
ROS1 6 117620577 117670578
6+
RET 10 43587031 43637031
7+
NTRK3 15 88458983 88508983
8+
ETV6 12 11997902 12047902
9+
EWSR1 22 29658122 29708122
10+
TMPRSS2 21 42841301 42891301
11+
FLI1 11 128652074 128702075
12+
HOOK3 8 42798356 42848357
13+
BRD4 19 15339962 15389963
14+
ATF1 12 51183062 51233063
15+
ETV1 7 13950462 14000463
16+
NUTM1 15 34615164 34665165

0 commit comments

Comments
 (0)