Rnaseq: GRCh38 iGenome

[higsch]

Is there a reason why GRCh38 is not included in the iGenomes conf file?

[apeltzer]

I think not - I'm actually finishing a grant application this week but will hopefully find some time next week to help out on adding missing stuff to iGenomes pretty soon 👍

Grch38 is there already, but I'm not sure whether the hisat2/star indices are!

[MaxUlysse]

Also, it's GRCh38 from NCBI, not from Ensembl

[marchoeppner]

You might also want to have a close look at http://lh3.github.io/2017/11/13/which-human-reference-genome-to-use before making a choice on this. Seems like the commonly used GRCh38 builds are generally not recommended for alignment-based analysis.

[ewels]

@MaxUlysse is currently working on adding the ensembl GRCh38 to AWS-iGenomes

[higsch]

On Uppmax, we already have it, right?

[ewels]

Related to https://github.com/ewels/AWS-iGenomes/issues/3

Basically, once we have ensembl GRCh38 available, we will update the config file to have this.

[adairama]

Dear all, may I ask when the ensembl GRCh38 is likely to be available please? Or if there is a workaround to this. Thank you.

[higsch]

You can use the --fasta and --gtf flag to use your own files. You can download them from http://www.ensembl.org/index.html.

If you then also add --saveReference the indices for STAR and the BED file is stored to results. If you want to rerun the pipeline in the same or another project you can add these generated references by --star_index and --bed12.

[adairama]

Thank you for this advice and this might be the best approach for now. Does anyone know which versions of Ensembl genome (non-masked primary assembly?) and gtf that igenomes used for GRCh37?

I prefer to use the same versions for GRCh38 which means I can minimize the need to re-run the pipeline when Ensembl GRCh38 is officially released on igenomes. I am not able to find this information from the igenomes website, their dumps or forums. I have emailed Illumina Technical support. I would be grateful for any pointers on this. Thank you.

[ewels]

Sadly I've never been able to find any further information from illumina about iGenomes beyond what is on the FTP (despite asking quite a few illumina employees who have asked around on my behalf).

See the main docs and the GRCh37 Annotation docs. The latter is pasted here for convenience:

The contents of the annotation directories were downloaded from Ensembl on: July 17, 2015.

Gene annotation files were downloaded from Ensembl release 75.
SmallRNA annotation files were downloaded from miRBase release 21.

[adairama]

Dear Phil, thank you for the links. It is a little sad/worrying that we are not entirely sure how such a widely used resource has been constructed. If they documented the process better, then the community could help contribute. I received a response from the Illumina Technical Support today. They are unsure of the timeline for Ensembl version of GRCH38. Further they said that

For the build method for Ensembl GRCh37, the information would be provided by the consortium during assembly (iGenome only curate the genome). You may want to check the information from Ensembl. For example, the GRCh38 from NCBI is indicated as haploid-with-alt-loci (https://www.ncbi.nlm.nih.gov/assembly?term=GRCh38&cmd=DetailsSearch )

Going by this logic, I arrived at this page: http://asia.ensembl.org/Homo_sapiens/Info/Annotation

I also had another look at the igenomes website and I found the following in one of their README

Sequence/Chromosomes/ contains individual chromosome sequences in fasta format.
Chromosomes are named according to their download source. Genome sequences from
UCSC are repeat-masked with lower-case characters; those from NCBI are not
repeat-masked; and those from Ensembl are repeat-masked with 'N' characters.
Masking was performed by the download source prior to download for the iGenomes.

Sequence/WholeGenomeFasta/ contains a multi-fasta files that contains all
chromosome sequences. This directory also contains samtools index file and
a picard sequence dictionary.

I am a bit lost as to whether these are useful information and how to reconcile them into something more meaningful.

[adairama]

BTW, the reason I am asking all of this is because the nextflow rnaseq pipeline with NCBI GRCh38 version fails for me as the corresponding GTF does not have biotypes column.

I could comment out lines 964 - 966 in main.nf or set the flag --fcGroupFeaturesType gene_name but sometimes I would like to see where the reads are aligning to.

My current workaround is to patch the NCBI GRCh38 GTF with biotype information from Ensembl GRCh38 GTF (downloaded from web) but this is a little ugly and I have gene name matching issues between the different version.

[higsch]

I used these ones:

Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
Homo_sapiens.GRCh38.92.gtf.gz

And it perfectly worked for me.

[drpatelh]

That post suggests that iGenomes Ensembl genomes were hard-masked...that's not great if that's the case. Good thing I tend to use the UCSC builds although these don't come with biotype information I think! I agree that they could be more transparent in terms of how the data is produced...

In terms of being consistent across assemblies, I personally don't think it matters that much as long as you consistently use the same reference for your analysis. Seeing as it doesn't look like this will be resolved any time soon I would just plug in your own GRCh38 reference to the pipelines and stick with that...

For future reference the following links will download the soft-masked primary assembly from Ensembl (release 96 in this case) and the associated GTF file:

ftp://ftp.ensembl.org/pub/release-96/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz

ftp://ftp.ensembl.org/pub/release-96/gtf/homo_sapiens/Homo_sapiens.GRCh38.96.gtf.gz

[adairama]

OK, looks like I been overthinking this. I will use the soft masked fasta and gtf as @drpatelh suggested. Thank you all.

[drpatelh]

So with a little testing it does indeed appear that the UCSC builds are soft-masked whereas the Ensembl and NCBI arent. I just grepped for a, c, g, t in the relevant fasta files from the same human assembly.

I then used galculator to work out the nucleotide frequency in each assembly:

| assembly | sequences | bases | A | C | G | T | N |
|----------------|-----------------|-------------|-----------|-----------|-----------|-----------|-----------|
| UCSC hg19 | 25 | 3095693984 | 844868045 | 585017944 | 585360437 | 846097277 | 234350281 |
| NCBI build37.1 | 25 | 3095693984 | 844868045 | 585017944 | 585360436 | 846097277 | 234350282 |
| Ensembl GRCh37 | 25 | 3095693982 | 844166209 | 584393900 | 584725573 | 845388980 | 237019320 |

The Ensembl reference does have many more Ns indicating that its hard-masked...

[sidorov-si]

Actually, it turns out that the GRCh37 from iGenomes is not masked at all: I just counted all different symbols in the fasta files of GRCh37 from iGenoms and downloaded from Ensembl (primary assembly: unmasked, soft-masked, and hard-masked). According to the following table:

Reference : Number of 'N'
GRCh37 (iGenomes) : 237,019,316
GRCh37 (original, unmasked) : 237,019,519
GRCh37 (original, soft-masked) : 228,862,183
GRCh37 (original, hard-masked) : 1,707,170,898

the number of 'N' in GRCh37 from iGenomes closely matches the number of 'N' in the original unmasked version.

[adairama]

Your table shows unmasked genome have 237 million Ns, while soft masked have 228 million Ns.

Shouldn't the soft masked genome have more Ns on account of masking?

[sidorov-si]

@adairama Good point! Interestingly, the total number of 'N' and 'n' in the soft-masked version of GRCh37 is 228,862,183 + 8,157,336 = 237,019,519, which is exactly the number of 'N' in the unmasked version. This implies that the Ensembl team produces the soft-masked version not from the unmasked one directly, but from a kind of internal 'raw' assembly. But, unfortunately, I don't know how exactly they produce their assemblies, so I may be wrong.

According to the GRCh37 README (ftp://ftp.ensembl.org/pub/release-75/fasta/homo_sapiens/dna/README), I regarded Homo_sapiens.GRCh37.75.dna.primary_assembly.fa as unmasked, Homo_sapiens.GRCh37.75.dna_rm.primary_assembly.fa as hard-masked, and Homo_sapiens.GRCh37.75.dna_sm.primary_assembly.fa as soft-masked. I used fasta-get-markov -m 0 -full -counts genome.fa command from MeMe-4.11.2 to count all symbols in these FASTA files.

[drpatelh]

I am not entirely sure what the fix is here! However, I have added a warning if using the NCBI GRCh38 genome from igenomes.config here https://github.com/nf-core/rnaseq/pull/459/commits/1e752cbe210fd680f8c9613046762817784f4979. This also auto-activates --skip_biotype_qc so the pipeline doesn't fail. The warning also has a link to this issue for information regarding the download of the Ensembl version of the assembly. Since we may be moving to RefGenie in the future, hopefully, it will become easier to manage genomes from different sources.