Skip to content

Instantly share code, notes, and snippets.

@anamariaelek
Last active May 15, 2021 22:09
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save anamariaelek/ba07d0f72280f907d3caa874f2a5cf2e to your computer and use it in GitHub Desktop.
Save anamariaelek/ba07d0f72280f907d3caa874f2a5cf2e to your computer and use it in GitHub Desktop.
Add cell barcode from read name to CB tag in bam file, subset nucleosome free reads, and save fragments file
#!/bin/bash
#$ -N cb
#$ -cwd
#$ -q long-sl7,mem_512_12h,mem_512
#$ -l virtual_free=50G,h_rt=12:00:00,disk=50G
#$ -M anamaria.elek@crg.eu
#$ -o cb.out
#$ -pe smp 16
#$ -j y
BAM=$1
SAM=${BAM%.bam}.sam
BAMCB=${BAM%.bam}.CB.bam
BAMCBSORT=${BAM%.bam}.CB.sorted.bam
BAMNFR=${BAM%%.bam}.CB.NFR.bam
BAMNFRSORT=${BAM%%.bam}.CB.NFR.sorted.bam
FRAGMENTS=${BAMNFRSORT%.bam}.fragments
NCORES=$2
echo "Adding CB tag to" $BAM
samtools view -H -@ $NCORES $BAM > $SAM
samtools view -@ $NCORES $BAM | awk '/^@/ {print;next} {N=split($1,n,":");print $0 "\tCB:Z:" n[1]}' >> $SAM
samtools view -bS -@ $NCORES $SAM > $BAMCB
samtools sort -@ $NCORES $BAMCB ${BAMCBSORT%.bam}
samtools index $BAMCBSORT
rm $SAM $BAMCB
echo "Subsetting nucleosome free reads from" $BAMCBSORT
samtools view -H -@ $NCORES $BAMCBSORT > $SAM
samtools view $BAMCBSORT | awk -v LEN=147 '{if ($9 <= LEN && $9 >= -(LEN) && $9 != 0 || $1 ~ /^@/) print $0}' >> $SAM
samtools view -bS -@ $NCORES $SAM > $BAMNFR
samtools sort -@ $NCORES $BAMNFR ${BAMNFRSORT%.bam}
samtools index $BAMNFRSORT
rm $SAM $BAMNFR
echo "Creating fragments file from" $BAMNFRSORT
sinto fragments -b $BAMNFRSORT -f ${FRAGMENTS}.tmp -t "CB" -m 30 -p $NCORES --use_chrom "(?i)^(scaffold)|(chr)"
sort -k 1,1 -k2,2n ${FRAGMENTS}.tmp > ${FRAGMENTS}
rm ${FRAGMENTS}.tmp
bgzip ${FRAGMENTS}
tabix -p bed ${FRAGMENTS}.gz
echo "Done!"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment