FLASH-seq UMI protocol V.3

# 0. Variables
BASECALL_DIR="/path/to/flowcell/Data/Intensities/BaseCalls/"
OUTPUT_DIR="/path/to/output_folder/"
SAMPLESHEET="/path/to/Demultiplexing_SampleSheet.csv"

# 1. Bcl2fastq
ulimit -n 10000
cd /path/to/flowcell/
bcl2fastq --input-dir $BASECALL_DIR --output-dir $OUTPUT_DIR --sample-sheet $SAMPLESHEET --create-fastq-for-index-reads --no-lane-splitting

zcat Undetermined_S0_I1_001.fastq.gz | awk -F' 1:N:0:' 'NR%4==1{print $2}' | sort | uniq -c > left_index.txt
sort -k1,1 left_index.txt

for file in ./out/*R1*
do
zcat $file | wc -l 
done

# 0. Variables
OUTPUTREF="/path/to/STAR_indexed_genome/"
FASTA="GRCh38.primary_assembly.genome.fa"
GTF="gencode.v34.primary_assembly.annotation.gtf"

# 1. Genome indexing
# sjdbOverhang should be adapted based on the read length (read_length - 1)
mkdir $OUTPUTREF

STAR --runThreadN 15 --runMode genomeGenerate --genomeDir $OUTPUTREF --genomeFastaFiles $FASTA --sjdbGTFfile $GTF --sjdbOverhang 99

# 1. Extract UMI in read 1
umi_tools extract --bc-pattern="^(?P<discard_1>AAGCAGTGGTATCAACGCAGAGT|AGCAGTGGTATCAACGCAGAGT|GCAGTGGTATCAACGCAGAGT|CAGTGGTATCAACGCAGAGT|AGTGGTATCAACGCAGAGT|GTGGTATCAACGCAGAGT|TGGTATCAACGCAGAGT|GGTATCAACGCAGAGT|GTATCAACGCAGAGT|ATCAACGCAGAGT|CAACGCAGAGT|AACGCAGAGT|ACGCAGAGT|GCAGAGT|CAGAGT|AGAGT|GAGT|AGT|GT)(?P<umi_1>.{8})(?P<discard_2>CTAACGG)(?P<discard_3>G{0,4})" --stdin=sample.R1.fastq.gz --stdout=umi.UMIinR1.R1.fq --read2-in=sample.R2.fastq.gz --read2-out=umi.UMIinR1.R2.fq --extract-method=regex



# 2. Extract UMI in read 2
umi_tools extract --bc-pattern="^(?P<discard_1>GAGT|AGT|GT)(?P<umi_1>.{8})(?P<discard_2>CTAACGG)(?P<discard_3>G{0,4})" --stdin=sample.R2.fastq.gz --stdout=umi.UMIinR2.R2.fq --read2-in=sample.R1.fastq.gz --read2-out=umi.UMIinR2.R1.fq --extract-method=regex

# 3. In very rare cases (<0.0001%) can get the UMI in both R1 and R2. Find read IDs with a UMI in both R1 and R2
cat umi.UMIinR1.R1.fq | uniq | awk 'NR%4==1{print}' | sed 's/\@//g' > names.R1umi.txt
cat umi.UMIinR2.R1.fq | uniq | awk 'NR%4==1{print}' | sed 's/\@//g' > names.R2umi.txt
cat names.R1umi.txt | sed 's/_........ .*$//g' > names.R1umi.cleaned
cat names.R2umi.txt | sed 's/_........ .*$//g' > names.R2umi.cleaned
comm -12 <(sort names.R1umi.cleaned) <(sort names.R2umi.cleaned) > R1R2.toFilterOut

echo "===> Number of R1-R2 with both a UMI: $(wc -l R1R2.toFilterOut) <==="
echo "===> Number of R1 UMI before cleanup: $(wc -l names.R1umi.txt) <==="
echo "===> Number of R2 UMI before cleanup: $(wc -l names.R2umi.txt) <==="

# 4. Filter out double UMI reads from UMI reads
grep -f R1R2.toFilterOut names.R1umi.txt > R1.toFilterOut
grep -f R1R2.toFilterOut names.R2umi.txt > R2.toFilterOut
$BBMAP_filter in=umi.UMIinR1.R1.fq in2=umi.UMIinR1.R2.fq out=umi.UMIinR1.R1.tmp out2=umi.UMIinR1.R2.tmp names=R1.toFilterOut include=f overwrite=t
$BBMAP_filter in=umi.UMIinR2.R1.fq in2=umi.UMIinR2.R2.fq out=umi.UMIinR2.R1.tmp out2=umi.UMIinR2.R2.tmp names=R2.toFilterOut include=f overwrite=t

mv umi.UMIinR1.R1.tmp umi.UMIinR1.R1.fq
mv umi.UMIinR2.R1.tmp umi.UMIinR2.R1.fq
mv umi.UMIinR1.R2.tmp umi.UMIinR1.R2.fq
mv umi.UMIinR2.R2.tmp umi.UMIinR2.R2.fq

echo "===> Number of R1 UMI after cleanup: $(grep -c \@ umi.UMIinR1.R1.fq) <==="
echo "===> Number of R2 UMI after cleanup: $(grep -c \@ umi.UMIinR2.R1.fq) <==="

# 1. Get Internal reads by excluding the UMI reads 
cat umi.UMIinR1.R1.fq | uniq | awk 'NR%4==1{print}' | sed 's/\@//g' | sed 's/_........//g' > names.R1umi.txt
cat umi.UMIinR2.R1.fq | uniq | awk 'NR%4==1{print}' | sed 's/\@//g' | sed 's/_........//g' > names.R2umi.txt
sed -i 's/_........//g' R1.toFilterOut
sed -i 's/_........//g' R2.toFilterOut
cat names.R1umi.txt names.R2umi.txt R1.toFilterOut R2.toFilterOut > names.umi.txt

$BBMAP_filter in=sample.R1.fastq.gz in2=sample.R2.fastq.gz out=internal.R1.fq out2=internal.R2.fq names=names.umi.txt include=f overwrite=t

echo "===> Number of Internal Reads after cleanup: $(grep -c \@ internal.R1.fq) <==="

# 1. Combine
cat umi.UMIinR1.R1.fq umi.UMIinR2.R2.fq > umi.R1.fq
cat umi.UMIinR1.R2.fq umi.UMIinR2.R1.fq > umi.R2.fq

# 2. Final Clean-up
rm R1.toFilterOut R2.toFilterOut toFilterOut.txt
rm names.* umi.UMIinR*.R*.fq

# 1. Trim Reads
bbduk.sh -Xmx48g in1=FASTQ/umi.R1.fq in2=FASTQ/umi.R2.fq out1=FASTQ/umi.R1.2.fq out2=FASTQ/umi.R2.2.fq t=32 ktrim=l ref=adapters.fa k=23 mink=7 hdist=1 hdist2=0 minlength=29 tbo
bbduk.sh -Xmx48g in1=FASTQ/umi.R1.2.fq in2=FASTQ/umi.R2.2.fq out1=FASTQ/umi.R1.trim.fq out2=FASTQ/umi.R2.trim.fq t=32 ktrim=r ref=adapters.fa k=23 mink=7 hdist=1 hdist2=0 minlength=29 tbo

bbduk.sh -Xmx48g in1=FASTQ/internal.R1.fq in2=FASTQ/internal.R2.fq out1=FASTQ/internal.R1.2.fq out2=FASTQ/internal.R2.2.fq t=32 ktrim=l ref=adapters.fa k=23 mink=7 hdist=1 hdist2=0 minlength=29 tbo
bbduk.sh -Xmx48g in1=FASTQ/internal.R1.2.fq in2=FASTQ/internal.R2.2.fq out1=FASTQ/internal.R1.trim.fq out2=FASTQ/internal.R2.trim.fq t=32 ktrim=r ref=adapters.fa k=23 mink=7 hdist=1 hdist2=0 minlength=29 tbo

# 2. Rename
mv FASTQ/umi.R1.trim.fq FASTQ/umi.R1.fq
mv FASTQ/umi.R2.trim.fq FASTQ/umi.R2.fq
mv FASTQ/internal.R1.trim.fq FASTQ/internal.R1.fq
mv FASTQ/internal.R2.trim.fq FASTQ/internal.R2.fq

# 3. Clean-up
rm FASTQ/internal.R1.2.fq FASTQ/internal.R2.fq FASTQ/umi.R1.2.fq FASTQ/umi.R2.2.fq

# 0. Variables
GENOME="/path/to/STAR_indexed_genome/"
FASTQ_R1="/path/to/umi.R1.fq”
FASTQ_R2="/path/to/umi.R2.fq"
ID=”sample_id”

# 1. Mapping
STAR --runThreadN 10 --limitBAMsortRAM 20000000000 --genomeLoad LoadAndKeep --genomeDir $GENOME --readFilesIn "$FASTQ_R1" "$FASTQ_R2" --readFilesCommand cat --limitSjdbInsertNsj 2000000 --seedSearchStartLmax 30 --outFilterIntronMotifs RemoveNoncanonicalUnannotated --outSAMtype BAM SortedByCoordinate --outFileNamePrefix STAR/"$ID"_UMI_


# 2. SAM to sorted BAM
# -F 260 filters out unmapped and secondary alignments
samtools view -@ 5 -Sb -F 260 "$ID"_UMI_Aligned.sortedByCoord.out.bam > "$ID"_UMI_Aligned.sortedByCoord.filtered.bam
samtools index "$ID"_UMI_Aligned.sortedByCoord.filtered.bam

# 0. Variables
GENOME="/path/to/STAR_indexed_genome/"
FASTQ_R1="/path/to/internal.R1.fq”
FASTQ_R2="/path/to/internal.R2.fq"
ID=”sample_id”

# 1. Mapping
mkdir STAR
STAR --runThreadN 10 --limitBAMsortRAM 20000000000 --genomeLoad LoadAndKeep --genomeDir $GENOME --readFilesIn "$FASTQ_R1" "$FASTQ_R2" --readFilesCommand cat --limitSjdbInsertNsj 2000000 --outFilterIntronMotifs RemoveNoncanonicalUnannotated --outSAMtype BAM SortedByCoordinate --outFileNamePrefix STAR/"$ID"_INTERNAL_

# 2. SAM to sorted BAM
# -F 260 filters out unmapped and secondary alignments
samtools view -@ 5 -Sb -F 260 "$ID"_INTERNAL_Aligned.sortedByCoord.out.bam > "$ID"_INTERNAL_Aligned.sortedByCoord.filtered.bam
samtools index "$ID"_INTERNAL_Aligned.sortedByCoord.filtered.bam

mkdir FEATURECOUNTS
 
# 1. Assign UMI reads to features
featureCounts -T 1 -p -t exon -g gene_name -s 1 --fracOverlap 0.25 -a "$GTF" -R BAM -o FEATURECOUNTS/"$ID"_Aligned.txt STAR/"$ID"_UMI_Aligned.sortedByCoord.filtered.bam
	
# 2. Sort and Index Reads
samtools sort  -@ 10 FEATURECOUNTS/"$ID"_UMI_Aligned.sortedByCoord.filtered.bam.featureCounts.bam -o FEATURECOUNTS/"$ID"_UMI_Aligned.sortedByCoord.filtered.bam.featureCounts.sorted.bam
samtools index FEATURECOUNTS/"$ID"_UMI_Aligned.sortedByCoord.filtered.bam.featureCounts.sorted.bam

# 3. Deduplicate and Count UMI Reads
umi_tools count --per-gene --paired --gene-tag=XT --chimeric-pairs=discard --unpaired-reads=discard --assigned-status-tag=XS -I FEATURECOUNTS/"$ID"_UMI_Aligned.sortedByCoord.filtered.bam.featureCounts.sorted.bam -S FEATURECOUNTS/"$ID".umi.counts.tsv.gz

mkdir FEATURECOUNTS
featureCounts -T 1 -p -t exon -g gene_name --fracOverlap 0.25 -a "$GTF" -o FEATURECOUNTS/"$ID"_INTERNAL_featureCounts.txt STAR/"$ID"_INTERNAL_Aligned.sortedByCoord.filtered.bam