Skip to content

Instantly share code, notes, and snippets.

@ag1805x
Created July 24, 2019 11:06
Show Gist options
  • Save ag1805x/55ba0d88f317f63423a4e54adc46eb1f to your computer and use it in GitHub Desktop.
Save ag1805x/55ba0d88f317f63423a4e54adc46eb1f to your computer and use it in GitHub Desktop.
#!/usr/bin/env bash
# make_rRNA.sh
# Kamil Slowikowski
# December 12, 2014
#
# Modified: Arindam Ghosh (July 24, 2019 )
#
#
# Referenc Genome: GRCh38.p5 Ensembl release 84
#
#
# 1. Prepare chromosome sizes file from fasta sequence if needed.
#
# ftp://ftp.ensembl.org/pub/release-84/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa
# cut -f1,2 Homo_sapiens.GRCh38.dna.primary_assembly.fa.fai > sizes.genome
#
# 2. Make an interval_list file suitable for CollectRnaSeqMetrics.jar.
#
# Ensembl genes:
#
# ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz
#
#
#
# Picard Tools CollectRnaSeqMetrics.jar:
#
# https://broadinstitute.github.io/picard/command-line-overview.html#CollectRnaSeqMetrics
chrom_sizes=sizes.genome
# rRNA interval_list file -------------------------------------------------
# Genes from Ensembl.
genes=/home/dell/Documents/Arindam/Work/ReferenceGenome/Human_84/Annotation/gtf/Homo_sapiens.GRCh38.84.gtf
# Output file suitable for Picard CollectRnaSeqMetrics.jar.
rRNA=GRCh38.p5.rRNA.interval_list
# Sequence names and lengths. (Must be tab-delimited.)
perl -lane 'print "\@SQ\tSN:$F[0]\tLN:$F[1]\tAS:GRCh38"' $chrom_sizes | \
grep -v _ \
>> $rRNA
# Intervals for rRNA transcripts.
grep 'gene_biotype "rRNA"' $genes | \
awk '$3 == "gene"' | \
cut -f1,4,5,7,9 | \
perl -lane '
/gene_id "([^"]+)"/ or die "no gene_id on $.";
print join "\t", (@F[0,1,2,3], $1)
' | \
sort -k1V -k2n -k3n \
>> $rRNA
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment