Cell Ranger ARC, printed on 11/11/2024
Steps to create the pre-built Cell Ranger ARC reference packages from the downloads page.
# Genome metadata genome="GRCh38" version="2020-A" # Set up source and build directories build="${genome}-${version}-build" mkdir -p "$build" # Download source files if they do not exist in reference-sources/ folder source="${genome}-${version}-reference-sources" mkdir -p "$source" fasta_url="http://ftp.ensembl.org/pub/release-98/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz" fasta_in="${source}/Homo_sapiens.GRCh38.dna.primary_assembly.fa" gtf_url="http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_32/gencode.v32.primary_assembly.annotation.gtf.gz" gtf_in="${source}/gencode.v32.primary_assembly.annotation.gtf" motifs_url="https://jaspar.genereg.net/download/data/2018/CORE/JASPAR2018_CORE_vertebrates_non-redundant_pfms_jaspar.txt" motifs_in="${source}/JASPAR2018_CORE_non-redundant_pfms_jaspar.txt" if [ ! -f "$fasta_in" ]; then curl -sS "$fasta_url" | zcat > "$fasta_in" fi if [ ! -f "$gtf_in" ]; then curl -sS "$gtf_url" | zcat > "$gtf_in" fi if [ ! -f "$motifs_in" ]; then curl -sS "$motifs_url" > "$motifs_in" fi # Modify sequence headers in the Ensembl FASTA to match the file # "GRCh38.primary_assembly.genome.fa" from GENCODE. Unplaced and unlocalized # sequences such as "KI270728.1" have the same names in both versions. # # Input FASTA: # >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF # # Output FASTA: # >chr1 1 fasta_modified="$build/$(basename "$fasta_in").modified" # sed commands: # 1. Replace metadata after space with original contig name, as in GENCODE # 2. Add "chr" to names of autosomes and sex chromosomes # 3. Handle the mitochrondrial chromosome cat "$fasta_in" \ | sed -E 's/^>(\S+).*/>\1 \1/' \ | sed -E 's/^>([0-9]+|[XY]) />chr\1 /' \ | sed -E 's/^>MT />chrM /' \ > "$fasta_modified" # Remove version suffix from transcript, gene, and exon IDs in order to match # previous Cell Ranger reference packages # # Input GTF: # ... gene_id "ENSG00000223972.5"; ... # Output GTF: # ... gene_id "ENSG00000223972"; gene_version "5"; ... gtf_modified="$build/$(basename "$gtf_in").modified" # Pattern matches Ensembl gene, transcript, and exon IDs for human or mouse: ID="(ENS(MUS)?[GTE][0-9]+)\.([0-9]+)" cat "$gtf_in" \ | sed -E 's/gene_id "'"$ID"'";/gene_id "\1"; gene_version "\3";/' \ | sed -E 's/transcript_id "'"$ID"'";/transcript_id "\1"; transcript_version "\3";/' \ | sed -E 's/exon_id "'"$ID"'";/exon_id "\1"; exon_version "\3";/' \ > "$gtf_modified" # Define string patterns for GTF tags # NOTES: # - Since GENCODE release 31/M22 (Ensembl 97), the "lncRNA" and "antisense" # biotypes are part of a more generic "lncRNA" biotype. # - These filters are relevant only to GTF files from GENCODE. The GTFs from # Ensembl release 98 have the following differences: # - The names "gene_biotype" and "transcript_biotype" are used instead of # "gene_type" and "transcript_type". # - Readthrough transcripts are present but are not marked with the # "readthrough_transcript" tag. # - Only the X chromosome versions of genes in the pseudoautosomal regions # are present, so there is no "PAR" tag. BIOTYPE_PATTERN=\ "(protein_coding|lncRNA|\ IG_C_gene|IG_D_gene|IG_J_gene|IG_LV_gene|IG_V_gene|\ IG_V_pseudogene|IG_J_pseudogene|IG_C_pseudogene|\ TR_C_gene|TR_D_gene|TR_J_gene|TR_V_gene|\ TR_V_pseudogene|TR_J_pseudogene)" GENE_PATTERN="gene_type \"${BIOTYPE_PATTERN}\"" TX_PATTERN="transcript_type \"${BIOTYPE_PATTERN}\"" READTHROUGH_PATTERN="tag \"readthrough_transcript\"" PAR_PATTERN="tag \"PAR\"" # Construct the gene ID allowlist. We filter the list of all transcripts # based on these criteria: # - allowable gene_type (biotype) # - allowable transcript_type (biotype) # - no "PAR" tag (only present for Y chromosome PAR) # - no "readthrough_transcript" tag # We then collect the list of gene IDs that have at least one associated # transcript passing the filters. cat "$gtf_modified" \ | awk '$3 == "transcript"' \ | grep -E "$GENE_PATTERN" \ | grep -E "$TX_PATTERN" \ | grep -Ev "$READTHROUGH_PATTERN" \ | grep -Ev "$PAR_PATTERN" \ | sed -E 's/.*(gene_id "[^"]+").*/\1/' \ | sort \ | uniq \ > "${build}/gene_allowlist" # Filter the GTF file based on the gene allowlist gtf_filtered="${build}/$(basename "$gtf_in").filtered" # Copy header lines beginning with "#" grep -E "^#" "$gtf_modified" > "$gtf_filtered" # Filter to the gene allowlist grep -Ff "${build}/gene_allowlist" "$gtf_modified" \ >> "$gtf_filtered" # Change motif headers so the human-readable motif name precedes the motif # identifier. So ">MA0004.1 Arnt" -> ">Arnt_MA0004.1". motifs_modified="$build/$(basename "$motifs_in").modified" awk '{ if ( substr($1, 1, 1) == ">" ) { print ">" $2 "_" substr($1,2) } else { print } }' "$motifs_in" > "$motifs_modified" # Create a config file config_in="${build}/config" echo """{ organism: \"Homo_sapiens\" genome: [\""$genome"\"] input_fasta: [\""$fasta_modified"\"] input_gtf: [\""$gtf_filtered"\"] input_motifs: \""$motifs_modified"\" non_nuclear_contigs: [\"chrM\"] }""" > "$config_in" # Create reference package cellranger-arc mkref --ref-version="$version" \ --config="$config_in"
# Genome metadata genome="mm10" version="2020-A" # Set up source and build directories build="${genome}-${version}-build" mkdir -p "$build" # Download source files if they do not exist in reference_sources/ folder source="${genome}-${version}-reference-sources" mkdir -p "$source" fasta_url="http://ftp.ensembl.org/pub/release-98/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz" fasta_in="${source}/Mus_musculus.GRCm38.dna.primary_assembly.fa" gtf_url="http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M23/gencode.vM23.primary_assembly.annotation.gtf.gz" gtf_in="${source}/gencode.vM23.primary_assembly.annotation.gtf" motifs_url="https://jaspar.genereg.net/download/data/2018/CORE/JASPAR2018_CORE_vertebrates_non-redundant_pfms_jaspar.txt" motifs_in="${source}/JASPAR2018_CORE_vertebrates_non-redundant_pfms_jaspar.txt" if [ ! -f "$fasta_in" ]; then curl -sS "$fasta_url" | zcat > "$fasta_in" fi if [ ! -f "$gtf_in" ]; then curl -sS "$gtf_url" | zcat > "$gtf_in" fi if [ ! -f "$motifs_in" ]; then curl -sS "$motifs_url" > "$motifs_in" fi # Modify sequence headers in the Ensembl FASTA to match the file # "GRCm38.primary_assembly.genome.fa" from GENCODE. Unplaced and unlocalized # sequences such as "GL456210.1" have the same names in both versions. # # Input FASTA: # >1 dna:chromosome chromosome:GRCm38:1:1:195471971:1 REF # # Output FASTA: # >chr1 1 fasta_modified="$build/$(basename "$fasta_in").modified" # sed commands: # 1. Replace metadata after space with original contig name, as in GENCODE # 2. Add "chr" to names of autosomes and sex chromosomes # 3. Handle the mitochrondrial chromosome cat "$fasta_in" \ | sed -E 's/^>(\S+).*/>\1 \1/' \ | sed -E 's/^>([0-9]+|[XY]) />chr\1 /' \ | sed -E 's/^>MT />chrM /' \ > "$fasta_modified" # Remove version suffix from transcript, gene, and exon IDs in order to match # previous Cell Ranger reference packages # # Input GTF: # ... gene_id "ENSMUSG00000102693.1"; ... # Output GTF: # ... gene_id "ENSMUSG00000102693"; gene_version "1"; ... gtf_modified="$build/$(basename "$gtf_in").modified" # Pattern matches Ensembl gene, transcript, and exon IDs for human or mouse: ID="(ENS(MUS)?[GTE][0-9]+)\.([0-9]+)" cat "$gtf_in" \ | sed -E 's/gene_id "'"$ID"'";/gene_id "\1"; gene_version "\3";/' \ | sed -E 's/transcript_id "'"$ID"'";/transcript_id "\1"; transcript_version "\3";/' \ | sed -E 's/exon_id "'"$ID"'";/exon_id "\1"; exon_version "\3";/' \ > "$gtf_modified" # Define string patterns for GTF tags # NOTES: # - Since GENCODE release 31/M22 (Ensembl 97), the "lncRNA" and "antisense" # biotypes are part of a more generic "lncRNA" biotype. # - These filters are relevant only to GTF files from GENCODE. The GTFs from # Ensembl release 98 have the following differences: # - The names "gene_biotype" and "transcript_biotype" are used instead of # "gene_type" and "transcript_type". # - Readthrough transcripts are present but are not marked with the # "readthrough_transcript" tag. BIOTYPE_PATTERN=\ "(protein_coding|lncRNA|\ IG_C_gene|IG_D_gene|IG_J_gene|IG_LV_gene|IG_V_gene|\ IG_V_pseudogene|IG_J_pseudogene|IG_C_pseudogene|\ TR_C_gene|TR_D_gene|TR_J_gene|TR_V_gene|\ TR_V_pseudogene|TR_J_pseudogene)" GENE_PATTERN="gene_type \"${BIOTYPE_PATTERN}\"" TX_PATTERN="transcript_type \"${BIOTYPE_PATTERN}\"" READTHROUGH_PATTERN="tag \"readthrough_transcript\"" # Construct the gene ID allowlist. We filter the list of all transcripts # based on these criteria: # - allowable gene_type (biotype) # - allowable transcript_type (biotype) # - no "readthrough_transcript" tag # We then collect the list of gene IDs that have at least one associated # transcript passing the filters. cat "$gtf_modified" \ | awk '$3 == "transcript"' \ | grep -E "$GENE_PATTERN" \ | grep -E "$TX_PATTERN" \ | grep -Ev "$READTHROUGH_PATTERN" \ | sed -E 's/.*(gene_id "[^"]+").*/\1/' \ | sort \ | uniq \ > "${build}/gene_allowlist" # Filter the GTF file based on the gene allowlist gtf_filtered="${build}/$(basename "$gtf_in").filtered" # Copy header lines beginning with "#" grep -E "^#" "$gtf_modified" > "$gtf_filtered" # Filter to the gene allowlist grep -Ff "${build}/gene_allowlist" "$gtf_modified" \ >> "$gtf_filtered" # Change motif headers so the human-readable motif name precedes the motif # identifier. So ">MA0004.1 Arnt" -> ">Arnt_MA0004.1". motifs_modified="$build/$(basename "$motifs_in").modified" awk '{ if ( substr($1, 1, 1) == ">" ) { print ">" $2 "_" substr($1,2) } else { print } }' "$motifs_in" > "$motifs_modified" # Create a config file config_in="${build}/config" echo """{ organism: \"Mus_musculus\" genome: [\""$genome"\"] input_fasta: [\""$fasta_modified"\"] input_gtf: [\""$gtf_filtered"\"] input_motifs: \""$motifs_modified"\" non_nuclear_contigs: [\"chrM\"] }""" > "$config_in" # Create reference package cellranger-arc mkref --ref-version="$version" \ --config="$config_in"