text:
Using PLINK to find SNPs in LD (Linkage Disequilibrium) with another set of SNPs
I have a list of SNPs (index SNPs). As my tagging SNPs, I want to use SNP data from here. The data has the following column headers: Chromosome, Position, MarkerName, Effect_allele,   Non_Effect_allele, Beta, SE, Pvalue.
My ultimate aim is to identify SNPs in the tagging set that are in LD with each of index SNPs. I need to calculate r2 between index SNPs and tagging SNPs within 1MB window and get a table of index vs tagging SNPs (r2>0.8).
How can I use PLINK to carry out this work? I'm new to PLINK and LD clumping

You can't generate LD statistics for summary data from SNPs, as the LD calculation relies on the comparison of variants at the genotype level. Plink expects data to be in a genotype matrix with either each marker for each line (probably better in your case), or each individual for each line. Plink v1.9, in particular, will accept VCF-formatted input files, when populated with 'GT' information that describes the variants for each individual.
If you're okay with using the 1000 genome population for calculating LD, plink v1.9-compatible datasets should be accessible from
here.
If not, you'll need to discuss your project with IGAP to see if you can have access to genotype-level data for LD testing. That's unlikely, because there are personal identification issues associated with non-summary data.
There's comprehensive documentation about the LD command(s) on the Plink 1.9 website here.

Thanks for your explanation @gringer. I will now use a plink v1.9-compatible dataset from 1000 genome  and then use plink  --noweb
--bfile 1000genomefile --clump SNPfile --clump-field P --clump-p1 1.0001 --clump-p2 1.0001 --clump-r2 0.8 --clump-kb 1000 --out resultfile. But I'm unsure if line-by-line listing of SNPs on SNPfile will work and what the correct format should be. Any suggestions?

EDIT:  --clump-p1 0.0001 --clump-p2 0.0001


text:
Transferring genomic features on new coordinates
I have a eukaryotic genome for which an updated sequence for a chromosome  was recently obtained. I want to map RNAseq reads on the genome (and perform other downstream analyses) and would like to use the most up-to-date information possible (so the 'new' sequence of that chromosome).
However, before performing the read mapping, I'd like to update the annotation of my genome (GFF3) to match with the 'new' coordinates of this chromosome.
Basically, I would like to transfer information from the 'old-coordinates' GFF3 to obtain a 'new-coordinates' GFF3, and if possible conserve all the information/hierarchy of the file (gene, mRNA, exon, etc.).
Example:
chrXX   source  gene    222 5942    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    222 5942    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    222 5794    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    5889    5942    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 222 5794    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 5889    5942    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

... should be updated to ...
chrXX   source  gene    333 6053    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    333 6053    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    333 5905    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    6000    6053    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 333 5905    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 6000    6053    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

The approach I tried (but I am not sure it is the way to do it):

Extract sequence of all the features that are going to be updated to fasta (bedtools getfasta)
Map these sequences to the genome with the updated chromosome (gmap with --nosplicing since sequences we are mapping correspond to genomic regions).
Create the updated GFF3 file. When doing that, I apply one supplementary rule: if a feature mapped on another chromosome, but there was also another alignment on the updated chromosome, prioritize the one that correspond to the same (updated) chromosome.

What would be the proper way to do such a thing?  Any suggestion of methods/tools is welcome!

I'd suggest to paste the full code used for each step, so people can advise you better.

I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:
mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

Translate psl files to chains in the directory chain:
mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

Merge short chains into longer ones into the directory chainMerge:
mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

concat and sort the chains:
cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

Need info about chromosome sizes for netting:
twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

Netting: identify alignable regions from chains:
mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

Finally, select the right alignable regions using the nets, creating a "liftOver" file:
netChainSubset net/all.net all.chain ci2ToCi3.liftOver

Run liftOver:
CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed

Thanks for the quick answer, will test it and report back.

In addition to the answer by heathobrien you might take a look at RATT: Rapid Annotation Transfer Tool. The approach is essentially the same, only using more modern approaches to for instance whole genome alignments. It uses nucmer, which when using the latest version (4.0) allows you to do whole genome alignments using parallel computing.

Are there any resources that show how to do this? The documentation seems really non-user friendly.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Integration of different microarray dataset to run GSEA
I'm planning to run a GSE Analysis on some microarray datasets. Until now what I have done is to pre-processed them using frma algorithm.
The reason I chose frma instead of gcrma or rma is that I wanted to merge these different dataset into one.
Here is how I'm planing to proceed:
i.e
dataset1 : 3x Controls , 3x Treatment1
dataset2 : 2x Controls , 2x Treatment2
dataset3 : 2x Controls , 2x Treatment3

After merging
7x Controls , 3x Treatment1 , 2x Treatment2 , 2x Treatment3

Then I will execute the GSEA for controls , Treatmen1 , Treatment2 and Treatment3 separably.
So the question is if it's more reasonable to follow the path that I described or to run the GSEA on the separate datasets (for example for each control ) and then get the intersect of the gene sets of these three repetitions.

Is there any reason you don't want to or can't do both?

What have you done to control that each dataset is (presumably) done under different conditions and thus have batch effects?  Are the controls the same in each dataset? Could you point me what is the frma algorithm ? It is the first time I read about it.

@Mitra : No there is no such a reason. I'm just asking to see if something is more reasonable or makes more statistical sense than the other.

@Llopis : At the merging stage I will remove the batch effects with comBat method (I could also try others as well). Here is a documentation for frma

Thanks, take into account that this will remove too the effect of each treatment. as it cannot distinguish between the batch effect and the treatment effect.

It makes more sense to evaluate by separate the pathway or gene set you want and see if in the three datasets result in a coherent message than to merge these datasets, as you will mix batch effects and treatment effects.
When in several datasets one looks for the same statistic it is usually considered meta analysis, not integration.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
How to specify resources for cluster in snakemake
I am trying to run snakemake version 3.5.4 on my cluster, where jobs are partially defined by a common template build according manual and also defined though rule-specific parameters. But for some reason cluster configuration is not expanding wildcards (and command does). When I run
snakemake data/genome.fa.gz.bwt -p --jobs 10 --cluster-config cluster.json --cluster "bsub -J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} -o {cluster.output} -e {cluster.error}"

-J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} expands as expected. But -o {cluster.output} -e {cluster.error} is replaced by literally logs/cluster/{rule}.{wildcards}.out and logs/cluster/{rule}.{wildcards}.err, i.e. there is no expansion in config file (the config file is attached bellow).
I am puzzled since exactly this lines are in the documentation. What am I doing wrong? Why snakemake is not expanding wildcards in the config file?
Here is my cluster.json
{
    "__default__" :
    {
        "queue"    : "bgee",
        "one_host" : "\"span[hosts=1]\"",
        "name"      : "{rule}.{wildcards}",
        "output"    : "logs/cluster/{rule}.{wildcards}.out",
        "error"     : "logs/cluster/{rule}.{wildcards}.err"
    },
}

Here is relevant part of Snakefile :
rule index_reference :
        threads: 1
        resources: M=20000000
        input : "data/{reference}/genome.fa.gz"
        output : "data/{reference}/genome.fa.gz.bwt"
        shell : "scripts/index_genome.sh {input}"

This feature was added in version 3.6.0. Look at release notes.
{wildcards} normally expands to a something that looks vaguely like a dictionary. In other words, something like key1=value1,key2=value2. Whether you really want to use that as a log file name is up to you.
Let's use a very contrived example to demonstrate this. Here is the Snakefile:
rule all:
    input: "input1.human.txt"

rule blah:
    output: "{sample}.{reference}.txt"
    shell: "touch {output}"

I'm using your cluster.json, so with that let's use a contrived cluster submission program:
snakemake --cluster-config cluster.json --cluster "echo 'rule.wildcards: {rule}.{wildcards} cluster: {cluster} -o {cluster.output} -e {cluster.error} -n {cluster.name}' >> log && " --jobs 1

If you run that and look in log, you can see that {rule} is expanded as expected, {wildcards} is expanded to sample=input1,reference=human. The output and error parts from cluster.json are expanded correctly for me. Below are the contents of log:
rule.wildcards: blah.sample=input1,reference=human cluster: one_host="span[hosts=1]",output=logs/cluster/blah.sample=input1,reference=human.out,name=blah.sample=input1,reference=human,queue=bgee,error=logs/cluster/blah.sample=input1,reference=human.err -o logs/cluster/blah.sample=input1,reference=human.out -e logs/cluster/blah.sample=input1,reference=human.err -n blah.sample=input1,reference=human

I'm not a big fan of using {wildcards} in paths unless you're doing a parameter sweep, but to each their own. Note that the above results are from snakemake-4.3.1.

The solutions are either, upgrade snakemake or just take all wildcards out of cluster.json file and write them directly to snakemake execution command :
snakemake -p --jobs 10 --cluster-config cluster.json --cluster "bsub -J "{rule}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} -o {"logs/cluster/{rule}.{wildcards}.out" -e "logs/cluster/{rule}.{wildcards}.err""

However, I ended up instead of using of json config file in writing a simple bash wrapper snakemake_cluster.sh that remembers the execution command for me in more readable format :
#!/usr/bin/env bash

snakemake $@ -p --jobs 10 --cluster "bsub \
-J {rule} \
-q bgee \
-n {threads} \
-M {resources.M} \
-R \"span[hosts=1]\" \
-o logs/cluster/{rule}.{wildcards}.out \
-e logs/cluster/{rule}.{wildcards}.err"

Edit : All spaces after --cluster command got to be escaped, otherwise it would add arguments to snakemake not to bsub.
Edit 2 : $@ means all arguments, therefore it can be either targets, flags or combinations (targets got to be specified first). So, snakemake_cluster.sh my_file builds my_file, snakemake_cluster.sh --quiet builds default target quietly, and snakemake_cluster.sh my_file --quiet builds my_file quietly.

As of the most recent snakemake (v5.10.0) the recommended way of defining resources is through a profile.
It seems you are trying to submit a job on an LSF cluster system - there is a snakemake profile for LSF to save you having to do this work yourself.
It reduces your snakemake submission to something like
snakemake --profile "lsf"

Luckily I don't have to deal with lsf anymore, but thanks anyway, I will take a look at this feature.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
What algorithms exist for reverse protein folding?
Are there any algorithms already in existence for reverse protein folding/docking (e.g. taking a desired surface map/docking profile and generating DNA which will produce a protein to fold into it?)
Protein folding seems straightforward (though computationally intensive) and protein docking is likewise straightforward, but thus far I haven't come across any protein unfolding algorithms, though I might just be searching for the wrong keywords.
Any advice in this area?

This is the holy grial of bioinformatics! Predict the structure of a protein from the amino acid chain is hard and has many errors usually, but doing the other way round is almost impossible AFAIK. But if you found something let me know :)

@Llopis just starting on it to be honest.  I've been wanting to make a useably abstracted programming language which compiles into actual organisms for a long time and always looked at it as something best approached by cross referencing semantic-style matching algorithms vs DNA with tagged trait data by species (and where possible by gene,) but I recently decided to go get a biotech degree and it elucidated much of what I was missing - it seems protein unfolding is the only step actually missing to employ standard programming techniques to get it done (plus we already have DNA printers.)

Read this article and the articles that reference it. It is an utopia to program and compile organisms that way. As one of my professors said "the paper holds everything but then you need to check it". Go step by step if you want to compile organisms.

Oh I'm pretty sure I have the details at this point, the issue is the protein unfolding - everything beyond that isn't much different from logic used in systems programming.  I'll check out the paper though, thank you.

The researchers behind FoldIt have the occasional design challenges that use human game players to help design proteins to fit a structure.
However, even with a known protein sequence there's still a lot of labwork involved in trying to find a good working DNA sequence. Just generating the DNA sequence is a difficult task due to redundancy in the genetic code (i.e. 64 codons for 20 amino acids + stop signals). When that is combined with cell-specific tRNA frequencies that can influence the timing of amino acid incorporation (and therefore how a protein folds as it is created), creating a sequence that reliably fits a shape ends up being quite difficult.


text:
How can I find out the code DRIMSeq uses to make dmplotProportions objects?
I'm using DRIMSeq to do some analysis but I'd like to make my own figures from the results ideally using what they've done as a base.
I'm trying to find how they make their ggplot objects in the function plotProportoins but when I use findMethods('plotProportoins') it ends up with a method dmPlotProportoins.
It says in the manual that it's a ggplot object but I can't seem to find how they make it. Any ideas?

You can look at the source code of the package and see how they create the object.

What's the best way to do it? findMethods gives me the source code of the function but how can I get the code of the whole package? (sorry if the answer should be obvious)

In the tar ball you have the source code of the package.

I didn't know that was a thing. Brilliant, thank you!

I was curious how they have done it as well, so here expanded @Liopis comment :
wget http://bioconductor.org/packages/release/bioc/src/contrib/DRIMSeq_1.6.0.tar.gz
tar -zxf DRIMSeq_1.6.0.tar.gz
# if you need to find it, or use grep
# find DRIMSeq -name "dm_plotProportions.R"
less DRIMSeq/R/dm_plotProportions.R

This is planned to change and be able to look at the code thorough a web interface


text:
Estimate the length of poly-A tails from randomly-primed RNAseq data
So a poly-A tail is a long chain of adenine nucleotides that is added to a messenger RNA (mRNA) molecule during RNA processing to increase the stability of the molecule.
For my project, I would like to estimate the length of poly-A tails from randomly-primed RNAseq data.
Any ideas how this could be done? Are there algorithms/tools which currently do this?
Perhaps I would need to model this somehow? All suggestions appreciated

Why do you have randomly primed RNAseq data? Which organisms are you working with (or it is a mix or a simulated data) ? What have you looked up for your project ? Is there any database of poly-A tail lengths maybe?

@Llopis HeLa cells.  But I would like to do this for a few other organisms as well. "Is there any database of poly-A tail lengths maybe?" It would have to be a data-driven measurement

Just look for polyA tracts at the end of sequences, and count them if they're larger than ~18 bases. I've done this with MinION cDNA reads by mapping the polyA adapter sequence (with an elongated polyA sequence) to the reads using LAST, and working out the length of aligned sequence.
Unfortunately, because the cDNA adapter finishes with a TVN sequence, rather than TTT, binding of the adapter is not necessarily a good representation of the true polyA length. I expect that Illumina adapters will have similar problems, if not worse (due to the read length limits).
It might be possible to use direct-RNA sequencing to better categorise the polyA lengths, but I've yet to look at the public direct-RNA datasets. Just in case you're interested in looking yourself, here's a direct RNA dataset for NA12878 from the wonderful nanopore-WGS consortium:
https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md
Here's the consortium's graph from that summary document showing polyA length:

"> It might be possible to use direct-RNA sequencing to better categorise the polyA lengths, but I'm not familiar with any available public datasets for that." There's no way to infer this from my own dataset?

You've done direct-RNA sequencing on the MinION?

No, I think I'm misunderstanding the reference to "public datasets"

Direct RNA sequencing is a very new technology that can only be done on a MinION. There's not much data available for it yet.

Can’t be done. If you already sequenced then I’m afraid the money is wasted (unless of course the data is good for something else).
The standard Illumina basecaller doesn’t deal well with homopolymers and will consequently call incorrect poly(A) lengths.
Amongst others, Narry Kim’s lab developed a dedicated method — TAIL-seq — to solve this problem. It’s a protocol for library prep using pulldown of biotin tagged tails (rather than random priming), plus dedicated basecalling software that uses spike-ins for calibration.
Be warned though: the protocol is reportedly technically challenging, and in order to use it you need an Illumina machine that gives you access to the raw image files before basecalling. Up-to-date Illumina machines won’t give you these files, full stop. So you’ll need to find a machine that runs an outdated firmware.

(This is a copy of my answer from Reddit; I only now noticed this question here.)

Thanks. I'll probably remove the Reddit question soon, so the link will be broken.


text:
How to submit a canu job on LSF high-performance computing cluster farm?
I am currently running canu on an LSF Linux server using the following script called assemble.sh:
canu -p species -d output genomeSize=100m -pacbio-raw subreads.fastq

I then submit this script with the following in a computer node:
bsub -q basement -E 'test -e /nfs/users/nfs_c/cr' -R "select[mem>25000] rusage[mem=25000]" -M25000 -o canu.o -e canu.e -J canu 'bash assemble.sh' -n8 -R"span[hosts=1]"

On my laptop, using 8 cores and 16 GB memory, assemble.sh took about 48h. On the cluster, it's been running for 72h and it is still running. I wonder if I am allocating resources correctly or if I can do something better. On my laptop, I can see in my activity monitor that canu is using many cores at the same time. Do you know how I can check this in an LSF server?
This page talks about submitting jobs on LSF. In particular, it mentions this: gridOptions="your options list" but I don't know how to use this option.
Any help would be appreciated.
checking the queue with bjobs -u all shows there are 125259 jobs running. My own job is running, it's not pending. Does the number of other jobs affect the running of my own job? I will ask IT if there are preferred writing locations. I did not consider that could change the speed of the job. Do you know a command to check how many resources my job is using?

I haven't used canu for several years but looking back at my labbook I ran it with two additional flags:
-minThreads=8 -maxThreads=8
My memory on the subject is a bit hazy but  think I had to play around with these settings to get it to work.  I also set the flag:
-maxMemory=6

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
dividing genome into non-overlapping windows using R
I have nearly 10 million SNPs located on 10 chromosomes. I want to divide the genome into non-overlapping windows of 15, 20 and 30 kb. Here is part of my SNP table:
head (sap_ids)
       snp_id    chr  pos
    Chr01__15043   1 15043
    Chr01__15079   1 15079
    Chr01__15139   1 15139
    Chr01__15165   1 15165
    ...
    ...
    ...
    Chr17__214708424  17 214708424
    Chr17__214708451  17 214708451
    Chr17__214708484  17 214708484
    Chr17__214708508  17 214708508
    Chr17__214708574  17 214708574

I have been using this code, but it gives me wrong output because each chromosome starts from position 1
win_size<-c(15000,30000,50000,100000)
res<- cbind(snp_ids,
       data.frame(lapply(setNames(win_size, paste("window",win_size, sep = "_")),
       function(x)as.numeric(ceiling(snp_ids$pos/x)))))

For instance, as I you see in the example below in window 4 I get SNPs from chromosomes 1, 3 and 5:
snp_id     chrom poistion   window
Chr01__58332   1 58332            4
Chr01__58335   1 58335            4
Chr01__58341   1 58341            4
Chr01__58450   1 58450            4
Chr01__58471   1 58471            4
Chr01__58530   1 58530            4
Chr01__58542   1 58542            4
Chr01__58641   1 58641            4
Chr03__45457   3 45457            4
Chr03__45604   3 45604            4
Chr04__56873   4 56873            4
Chr04__57387   4 57387            4
Chr04__57399   4 57399            4
Chr04__57528   4 57528            4
Chr04__58419   4 58419            4
Chr04__59670   4 59670            4
Chr04__59704   4 59704            4

I guess the solution would be to run the code above by looping through each chromosome I tried something like this:
for (i in 1:10){
       res <- cbind(cats,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 paste(cats$chr[cats$chr == i], ceiling(cats$pos/w), sep = "_"))
               ))
       }

But still not working appropriately (column 4 should be 1_1, for the chr 1)!
Chr01__912    1  912          2_1
Chr01__944    1  944          2_1
Chr01__1107   1  1107         2_1
Chr01__1118   1  1118         2_1

Can anyone help me to figure out how can I modify the code to create these non-overlapping windows properly?

Why do you have 4 fields in sap_ids for chr17 and 3 for chr1? Why does that change? Do your positions increase with the chromosome numbers (are all positions on chr2 greater than those on chr1) or do you count each chromosome separately (so each chromosome starts from position 1)? If the latter, it makes sense that your approach fails.

I do not get the part with 15, 20 and 30 kb windows. You want the windows to be of various sizes dependent of number of SNPs inside or you want to have one window size that will correspond on average to ~200 SNPs?

@terdon, I just edited my post slightly. see my edits. yes, actually in my case I count each chromosome separately which means each chromosome starts from position1. Any idea how to do it in R?   @ Kamil S Jaron, I just want to do my analysis with different window sizes! that's it.

@Anna To avoid removing the answers I rolled back to the previous version. If your question is not fully answered please [edit] it or comment explaining why not. You can always deselect the accepted answer.

Please stop radically changing your question! If you need more help, ask a new question explaining what else you need. But since this has been answered, changing it renders the answers obsolete.

The question is a bit confusing to me, at it's core I understand you want non-overlapping windows across chromosomes.
One way to achieve is to use GenomicRanges::tileGenome function, which needs the chromosome lengths as input, e.g.
library(Biostrings)
mygenome <- readDNAStringSet(list.files(mypath,"mygenome.fa$",full=TRUE))
chrSizes <- width(mygenome)
names(chrSizes) <- names(mygenome)
print(chrSizes)
#     Chr1     Chr2     Chr3     Chr4     Chr5
# 18585056 19698289 23459830 26975502 30427671

This can than be fed into tileGenomes:
library(GenomicRanges)
bins   <- tileGenome(chrSizes, tilewidth=5e5, cut.last.tile.in.chrom=T)

This creates a GRange object which gives you a whole repertoire of post processing options, i.e. overlapping with SNPs or NGS reads (I use the latter for genome wide coverage plots):
print(bins)
# GRanges object with 240 ranges and 0 metadata columns:
#         seqnames               ranges strand
#            <Rle>            <IRanges>  <Rle>
#     [1]     Chr1   [      1,  500000]      *
#     [2]     Chr1   [ 500001, 1000000]      *
#     [3]     Chr1   [1000001, 1500000]      *
#     [4]     Chr1   [1500001, 2000000]      *
#     [5]     Chr1   [2000001, 2500000]      *
#     ...      ...                  ...    ...
#   [236]     Chr5 [28000001, 28500000]      *
#   [237]     Chr5 [28500001, 29000000]      *
#   [238]     Chr5 [29000001, 29500000]      *
#   [239]     Chr5 [29500001, 30000000]      *
#   [240]     Chr5 [30000001, 30427671]      *
#   -------
#   seqinfo: 5 sequences from an unspecified genome

Relatedly, there's the tile() command from GenomicRanges that will produce a GRangesList object of the input. It won't produce exactly what was desired by OP, but it may be a more convenient alternative.

Thanks so much @Sebastian Müller. yes. This is actually a better approach than writing code for that.

Glad it is useful! Maybe you could ask even ask your SNP question in a different question and take it out here so the the actual question becomes more in-line with the title (which is almost self-explanatory by itself)?

Also for completeness, there is also the GenomicRanges::slidingWindows function which takes a GRanges object designed to results in overlapping windows, but can be made to produce non-overlapping windows as well by setting width = step


text:
Automating sequence download from Wormbase
I'm working on the C. elegans model organism, for which there exists an online biological information database: Wormbase.
I'm able to manually search for a transcript (say "K06C4.12") and this gets me the following URL:
http://www.wormbase.org/species/c_elegans/cds/K06C4.12
Say I have this transcript ID: "F08G2.2". I guess that I can form a valid similar URL, as follows:
http://www.wormbase.org/species/c_elegans/cds/F08G2.2
Indeed, it works as expected.
Now, I would like to go further and be able to download in an automated way what can be obtained by clicking on one of the "download" buttons in the "Sequence" section of the above pages. More specifically, I'm interested in the "spliced + UTR" and "conceptual translation" fasta sequences.
If I use the DOM inspector in Firefox (Shift-Ctrl-C), I can locate things like the following:

    download

        >K06C4.12 spliced + UTR
        ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTCACCAAGCCAAAGGACGGAAAGAAGAG
        ACGTCATGCCCGTAAAGAATCCTACTCCGTGTACATCTACCGTGTGCTCAAGCAGGTTCATCCAGACACTGGAGTCTCTT
        CCAAGGCTATGTCCATCATGAACTCTTTCGTCAACGACGTCTTCGAGCGTATCGCTGCTGAGGCATCTCGTCTTGCTCAC
        TACAACAAGAGATCCACCATCTCTTCTCGCGAAATTCAAACTGCCGTCCGTTTGATTCTTCCAGGAGAGCTTGCCAAGCA
        CGCCGTCTCTGAGGGAACCAAGGCTGTCACCAAGTACACCTCTAGCAAGTAA

K06C4.12conceptual translation.fasta

        >conceptual translation
        MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH
        YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

    view conceptual translation (123 aa)

    >conceptual translation

    MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH

    YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

I don't find these if I simply download the page with wget and search in the source code, so they must be dynamically generated.
How could I proceed to get the sequences?
The documentation about Wormbases's REST API is not very detailed, so I don't understand how to form an URL that would give me the desired sequences.

Did you try to download the whole database and search there? (I am not sure if you need large amounts of data or not ...)

I don't know where to download it, neither how it is structured. I don't need a large amount of data, so I hope querying the site should be OK.

To download it they said it is here. There should be a ftp site, so you could try that too.

Indeed, they even provide a scary shell script (first time I see a "goto" instruction in a shell script) to install the database: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/acedb/INSTALL
I'm not sure I'm ready to delve into this for the moment...

Further inspecting the source code of the downloaded page, I found the following:

    <span class="ui-icon ui-icon-close module-close" wname="sequences"

    Sequences

This suggested me to have a look at http://www.wormbase.org/rest/widget/transcript/K06C4.12/sequences and I obtained a page containing the elements found in the DOM inspector.
This can be parsed, for instance using the bs4 python module:
# To get the page
import requests
# To parse the page content
from bs4 import BeautifulSoup
# To get this into Biopython
from io import StringIO
from Bio import SeqIO

rest_template = "http://www.wormbase.org/rest/widget/transcript/{accession}/sequences"
rest = rest_template.format(accession="K06C4.12")
# This did not work with "html.parser"
soup = BeautifulSoup(requests.get(rest).text, "html5lib")

def generate_fasta(soup):
    name = None
    for div in soup.find_all("div"):
        if "ui-helper-hidden" in div["class"]:
            if div["id"] == "content":
                fasta = SeqIO.read(StringIO(div.text), "fasta")
                description_fields = fasta.description.split()
                if len(description_fields) == 4:
                    fasta.id = "_".join(description_fields[:2])
                    if name is None:
                        name = fasta.name
                elif len(description_fields) == 2:
                    # This assumes that the "conceptual_translation"
                    # is after the "unspliced + UTR" or "spliced + UTR"
                    # Otherwise, one would need to pass
                    # the accession as function parameter
                    fasta.id = f"{name}_translation"
                else:
                    fasta.id = "_".join(description_fields)
                yield (fasta.id, fasta)

records_dict = dict(generate_fasta(soup))
print(records_dict['K06C4.12_spliced'])
print(records_dict['K06C4.12_translation'])

Output:
ID: K06C4.12_spliced
Name: K06C4.12
Description: K06C4.12 spliced + UTR
Number of features: 0
Seq('ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTC...TAA', SingleLetterAlphabet())
ID: K06C4.12_translation
Name: conceptual
Description: conceptual translation
Number of features: 0
Seq('MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSS...SK*', SingleLetterAlphabet())


text:
Identifying relevant SNPs from a list
I have a list of all SNPs (a list of RsIDs) that fall in the coding sequence of several thousand human genes, and I'm looking for a convenient method to come up with a list of known SNPs that alter an amino acid in those genes.
I'll be glad to hear your thoughts.

Are you open to commercial (non-free) solutions? Also, http://varso.me is a great tool but it can't take batch requests.

What you are looking for is SNP annotation. If you have the chromosome:position reference and alternate alleles for your SNPs of interest, it can be as simple as uploading them to the variant effect predictor.
http://grch37.ensembl.org/Homo_sapiens/Tools/VEP
This will give you the predicted protein change and novelty of the variant with respect to known variation in large scale sequencing projects like the 1000 genomes.
If you cannot afford to upload data on the web, there are downloadable command line software like SNPeff, Variant effect predictor etc but it would require some more bioinformatics knowledge of converting your data to a VCF format and obtaining a human reference genome.

ANNOVAR is another tool that will help you functionally annotate genetic variants. It will tell you if the variant is known in COSMIC or Clinvar databases, which AA change will occur, if it's a frameshift, synonymous, non-synonymous, etc...


text:
Joining scaffolds to get full genome?
I am using Ray to assemble a about 7.2 kb picornavirus genome. I am using k-mer of 55. Following is the exact code that I used.
mpiexec -n 1 Ray -k 55 -s 2-CA-SVA-2_S2_L001_R1_001.fastq_q30_trim_LR_15.fasta  -o rayd55tr

This worked quite well as I got a scaffolds.fasta file with two scaffolds one with 5730 bp and another with 1548 bp.
Now, I am confused how I can join those scaffolds to get one linear genome. I joined them manually and blasted just to check and I got something like below :.
Is there any tool to join these two contigs together to have one complete linear genome? I have reference genome too if needed.

What did you blast against? Your image just shows a few identical sequences, but we can't interpret any blast results without knowing what was blasted against what. Also, if you have a reference genome, why didn't you use that to assemble the contigs?

Are there any good tools available for reference based contigs joining ?

It is quite important to know what kind of sequencing data you have used for people to provide accurate tips.

Ragout worked well for me. It is a reference-based scaffolder.
ragout.py --outdir output/ input/ragout.rcp > message.txt 2> error.txt --threads 8 --overwrite

The ragout.rcp file contains a 'recipe' with the reference sequence and the contigs in FASTA format.:
.references = ensembl
.target = assembly

ensembl.fasta = ref.fa
assembly.fasta = contigs.fa

In my case, it reduced the number of contigs 10-fold and increased the N50 10-fold. It added N's between the contigs.

Please consider expanding your answer. How/when does it work well ? How will it help the OP in their case?

I have added my code and what Ragout did for me. What does 'OP' mean?

Great [edit]! OP comes from Original Poster to refer the one that asked the question. I hope this will be useful.

In some cases you might consider reference-based scaffolding, as per charlesdarwin's answer. I would argue that in most cases this will only provide limited information, as this will only really provide you with the order of the contigs as it is in your reference, as opposed to how it is in your sampled organism. In that sense it really only provides pseudo-scaffolding that makes your assembly look artificially good.
In all cases you will have to check wether your raw data supports your scaffolding solution.
Example: you will never be able to call structural variants across contig boundaries if you use reference based scaffolding.
I would suggest you try different k-mer sized for your assembly algorithm, and look into dedicated scaffolding tools. Without knowing your type of sequencing data it is hard to suggest anything specific.


text:
Estimating BAM file from compressed fastq file size
Is there a way to estimate the size of a BAM file will have after mapping with BWA?
The input file are two mates fastq files, compressed with gzip, each one about 70G.

A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference.
If you are concerned about space and have a reference sequence available, you can store mapped reads as CRAM format, which will save a lot more space.

Could also be smaller if you remove duplicate reads. I doubt it will ever be significantly different in size from the original fastq file, but a +/- a few hundred megabytes (more for larger files) or so is to be expected.

A sorted BAM will be (potentially substantially) smaller than an equivalent FASTQ file, even when containing mapping information.


text:
How can I download genePred file from UCSC?
I need a gtf annotation file. UCSC doesn't give us a proper gtf file with distinct gene_id and transcript_id. It asked us to get a genePred file to convert to gtf.
But where can we get genePred files?

It turned out that using table browser on UCSC we could download an output with all fields. That output actually contains a GenePred file but has one extra column on the left end.
The genePred file format definition is here.


text:
Gene function annotation - bacterial genome
I'm trying to annotate a genome to find all genes with a specific function. I have a FASTA and the read FASTQs - I'd like to assign the functional group of the identified proteins (e.g. Kegg orthology) automatically.
For more context, I have whole genome sequences of bacteria - I'm trying to find all genes involved in membrane transport. I'd either like a tool to retrieve gene ontologies from a list of gene names (input: genA output: Membrane Transport), or a tool that outputs the function as part of the annotation pipeline (e.g. genbank output with GO function defined).
Or posed in a different way, how would you identify all genes with a specific function from prokaryotic whole genome sequence data?
Are there any scripts that can do this?

Please [edit] your question and give us more details about what you have and what you need. Do you have gene models for the bacterium? I mean, do you first have to find the genes, then the gene functions, or do you already have the genes identified? You have FASTA and FASTQ of what? The entire genome? Raw reads? Are they assembled? Is there a reference sequence for the bacterium?

Thanks for your reply - edited in a few more details. I have the raw read data, and relative de novo assemblies of whole genome prokaryotes. I have a few species I'm looking at, a few E. coli, Staph, and Yersinia. Ref seqs available for all, but I'll be annotating draft clinical genomes.
I'm looking for a way to find all genes present in an assembly that match a specific function (e.g. membrane transport). In the genbank format for example I see "/function=" which is exactly what I want to assign - but I can't find a tool that writes that data. Prokka for e.g. does not write this by default

The PATRICBrc.org is a free resource you can annotate your genome and compare it with publicly available genomes.
Refer this tutorial for more detail.
Once you have annotated your genome in PATRIC, you will have PATRIC Local families (k-mer based protein family), which will help you identify all genes with a specific function.

Im not 100% on what you are looking for - you seem to jump from annotation to detection of genes in your question.
But if you have a database of genes you want to detect (presence/ absence) you could try srst2
which basically takes your raw reads (fastq) and maps them to the database.  it will give you a report on any hits, and any snps/ indels in the hit. The only problem is making you own database is difficult as you must adhere to the fasta header format.


text:
Scaffolding a genome with hybrid data
I am assembling a ~500MB genome, and have ~150x long reads and ~200x 150bp PE short reads, with a ~400bp insert size.
I've done a lot of work with minimap+miniasm, and have what I think is a good set of unitigs and a .gfa file from miniasm. Depending on the parameters, I have a few hundred to ~1K unitigs. I can easily polish these unitigs with pilon, racon, and nanopolish, so that they are as error free as possible.
The question I have is: given access to

Polished (or not) unitigs
150x long reads (long enough to span almost all repeats in the genome)
200x 150bp PE short reads (~400bp insert size)

What would you do to scaffold a genome? And what are the pros and cons of your suggested approach?
What I need here is a draft genome. I would prefer a conservative approach that is more fragmented and has fewer misassemblies. So of course one option is just to keep the unitigs and make no attempt at scaffolding. But it does strike me there may be suitably conservative approaches I haven't considered. E.g. perhaps I could use ABYSS to scaffold my polished unitigs using my short-read data, thereby leveraging whatever information I can from the PE short reads?

You don't specify how long the long reads are, but if they are longer than 700 bp, then the PE reads won't be much help for scaffolding as they will contain less long-range info than the long reads. The short reads will probably be most useful for polishing the long reads.

Good point. The long reads are almost all >10KB. I guess I had a vague notion that there might be some useful information in the PE reads.

Try LR_Gapcloser.
I've used L_RNA_Scaffolder for trying to scaffold a genome (which turned out to be more complex than I had expected). It looks like LR_Gapcloser (written by the same people) is similar, but designed specifically for scaffolding using long-reads.
That page also suggests PBJelly and GMCloser as competing tools.


text:
How to get Nanopore MinION fast5 from SRA
I found some Nanopore MinION data on SRA, which I would like to investigate. I use sratoolkit for Illumina data all the time, but I am not sure how to get the fast5 file from the sra file.
What I have tried is to prefecth it, but then I don't know which command to use for .sra to .fast5 conversion. Not with fastq-dump, I assume?
prefetch -v SRR5286963 --max-size 50000000

When I was trying to upload MinION data to SRA, they wouldn't allow me to upload only the raw signal data (i.e. FAST5 files). Every file needed to contain a FASTQ sequence, and this sequence was extracted from the file and stored in the archive. Given this, I don't expect it will be possible to retrieve signal or event information from SRA. I also don't expect this will be possible in the future either, given that the SRA database expects every record to have a FASTQ-like sequence attached.
However, EBI/ENA do allow the upload of raw FAST5 files to the database (e.g. here), but there are no FAST5 reads for the specific single-cell B1a study you have indicated. If you're interested specifically in that study, then I recommend that you ask Miten Jain directly for access to those reads; he should have another method to give you access to them.
ENA is "ready" for ONT reads, but researchers need to upload their reads to ENA in order for that database to be useful.

Thanks for sharing that, I guess SRA and ENA are not completely ready yet for ONT.

I would suggest trying an alternative way, as the FAST5 is a fairly new entry in SRA - 6 public datasets as of today.
You can search the ERX/ERR entry you got from SRA in ENA and either (1) direct download it via ftp or (2) redirect to Galaxy.
Here's a tutorial that might be interesting for you: basically, you first get the archive from ENA, then find the FAST5 files in the top level of the downloaded archive, all to be processed in R using Fast5Summary objects from IONiseR package.

Thanks, great idea! However, on ENA I can only download fastq or the SRA file or am I missing something?

I guess there are no fast5 associated with SRR5286963 (simple filtered search), so maybe the easiest solution would be to directly ask the publishers - there is possibility that they just uploaded the fastq (e.g. from poretools extraction) instead of the fast5, because of the archive size.


text:
differential expression multiple studies
I was wondering if we can perform differential expression on samples from different library and read lengths.
What we intend to do is,

n = 6 tumors (recently sequenced using Hiseq 2500 read length 150, paired end)
n = 5 normals (previously sequenced using Hiseq read length 100, paired end, these are not matched normal)
n = 100 ( normals, from tcga read length 100, paired)

What we would like to do is compare n=6 Vs 100 + 5 normals. What steps should we be considering before we move ahead with this analysis?

It might be nice to add what steps you have considered already. For instance: how will you be able to tell if the observed differences you find are due to cancer/normal instead of due to different sequencing machines, different labs that produced the data or different people extracting the RNA?

I have no clue we still at design stage, what you have asked is a good point that we need to consider and look an answer for.

Adding experiments 1 & 2 should be fine without modification on the computer side of things, as long as your sample preparation is consistent. To be ultra-cautious, you could add in the experiment as an additional covariate for the differential expression model.
I'm less confident about adding in TCGA data; it's probably going to look better than your own samples, which might cause problems in calculating differential expression from raw reads. You can do a separate normalisation/correction of the TCGA data, but that will limit the inferences that can be made about differential expression -- without raw count information, the shot noise associated with read sampling cannot be included in the statistical model. I notice that TCGA has already done their own filtering / expression, and that would be preferable to use (mRNA sequencing, data subtype 'Gene'):
https://tcga-data.nci.nih.gov/docs/publications/tcga/datatype.html
However, even with data that has been properly normalised and filtered to exclude batch effects, there will still be population differences that influence results. My recommendation would be to do some subsampling to reduce the population effects, something like at least the following comparisons:

6 Vs 100 + 5 normals
6 Vs 100 normals
6 Vs 50 + 3 normals [randomly sampled from the two normal groups]

And only consider differential expression that is consistently observed in all three of these comparisons. Note that the third "subsampled" dataset can be generated more than once, taking different subgroups each time.
If you are going to use separately normalised data, then DESeq2 may not be appropiate (because it expects count-level data). Something like Limma could work for DE analysis, see chapter 15 of the user's guide for RNASeq analysis. Unfortunately I'm not very familiar with using Limma, so can't provide any advice on things like what to watch out for in the input data and/or results.

since I have not analyzed such data what sort of normalization would you reckon, before we push this to DESeq2

any suggestions how can i compare this?


text:
Mirbase and blast locally
I am trying to reproduce locally on my computer what I get running mirbase on their website using BLAST. The search sequences is: mature miRNAs which I had downloaded on my computer and make it as a BLAST database with command:
./makeblastdb -in /home/marianoavino/Downloads/mature.fa -dbtype 'nucl' -out /home/marianoavino/Downloads/mature

then on mirbase I see they use an e-value of 10, which I leave locally.
On mirbase at the end of the analysis they give you these parameter setting:
Search parameters

Search algorithm:   BLASTN
Sequence database: mature
Evalue cutoff: 10
Max alignments: 100
Word size: 4
Match score: +5
Mismatch penalty: -4

and this is the command line I use on my computer for BLAST
./blastn -db /home/marianoavino/Downloads/mature -evalue 10 -word_size 4 -query /home/marianoavino/Downloads/testinputblast.fasta -task "blastn" -out /home/marianoavino/Downloads/testBLast.out

The results of the two analysis are different, with mirbase finding much more stuff than local BLAST.
Do you have any idea on which parameters I should use on local blast command line to match those listed mirbase parameters in order to get the same answer?

You don't mention which version of BLAST you have installed locally and if that matches the web version.  Can you also provide a link the web version of mirbase?

BLAST version on my computer is 2.7.1+, the link for mirbase is http://mirbase.org/search.shtml, not sure which version of Blast they use. Thanks.

I would start with finding out what version of BLAST they use and then make sure you have the same one.  Then I would look at the man pages for BLAST and get the match score and mismatch penalty set.  Also even though it doesn't seem like it should matter the max alignments setting can change the results so I would have BLASt do that filtering for you as well.


text:
How can I create my own GO slim?
I know there are existing GO slims predefined. But does anyone know how I can create my own collection of GO terms for enrichment analysis?

I have never created one, but a GO slims is just a list of GO associated with genes, so you just need to download the list of GO and associate them with your list of genes.  What have you tried/searched? What do you need that it is not covered in the FAQ?

GO slims are just a subset of the gene ontology graph. Just a group of terms, no genes. We can't help you create one if you don't explain what terms you want and what sort of analysis you want to do.

I have a list of genes enriched in kidney. So I want to get a list of kidney-related terms and see which of them are significantly enriched.

@Dr_Hope Then what you really want is to test an enrichment in GO related to kidney  not create a GO slim. But if you only use a subset of genes and a subset of GOs the test would be biased (AFAIk) if you don't tell the suitable background population of genes and GO terms.

I've already done a canonical GO to check the background. The enriched terms are mostly kidney-related because our sample is a kidney. That's why I need a GO slim to check kidney-only terms.

authors: กิ๊ฟเลิกล่ะสายแข็ง


text:
Is there a more elegant solution to the bwa-mem: paired reads have different names error?
I'm currently trying to run bwa-mem on Influenza substrains using the following command:
~/bwa mem h5n1_1_cons.fa h5n1_1_read1.fq h5n1_1_read2.fq

h5n1_1_cons.fa is the consensus sequence for substrain h5n1_1, and the fq files are paired end reads 1 and 2. Running the program gives me this error:
[mem_sam_pe] paired reads have different names: "JQ714243-920-1", "JQ714243-920-2"

Visual inspection shows that the read names are identical until the last character. I have found a package (bbmap) with a utility, bbrename.sh, that renames the sequences, producing:
./bbrename.sh in=h5n1_1_read1.fq out=fixed prefix=h5n1_1

which yields the following output
-bash-4.2$ head fixh5n1_1_read1.fq
@h5n1_0
gaaaatgtctcaggagttgatgcatggtgtttaatctttggtttgcgctg
+
@4>C5?7FA?9@@@;1,-4B+E??4E9,DA>6CB)*2.)?7/.'65+75/
@h5n1_1
catcagttggattttgcctaaggatgtccaccatccttatcccaccaatt
+
7(9?)=C3AEF+@@F=C-E>EE57CD:;'?.*?7;>:3>2.;29..)=2'
@h5n1_2
catgcatttgttggaataatgttctcacaaatcaactgtattgaccgctg

This fixes the problem and allows BWA-MEM to work. However, header information is discarded. Is there another way to ensure the header names are the same while retaining all header information?

Please [edit] your question and show us a few examples of the read names in your original input files. The output of awk 'NR % 4 == 1{print $1}' h5n1_1_read1.fq | head and awk 'NR % 4 == 1{print $1}' h5n1_1_read2.fq | head for example.

Since you mention bbmap already, why not use that for mapping your reads altogether? It is after all called bbmap

Assuming the read name lines only contain things like @JQ714243-920-1 and nothing else, then sed -i '1~4 s/-[12]$//' file.fastq (N.B., example updated a bit for safety, thanks to @terdon) will remove the -1 and -2 suffixes.
However, you should try to figure out why the read names were munged into that form to begin with (at least Illumina machines won't produce that, I don't know about others).

You could also replace the the second - with a space. Most mapping tools only regard the part of the id before the first space as the read name.


text:
Duplicate genes with RSEM counts: Which one to choose?
I have Ensembl ids in the first column and samples with RSEM counts data in other columns. I converted Ensembl ids to gene symbols. Now I see there are three genes repeated twice.
Gene   S1       S2       S3
COG8  804.07  1475.16  323.80
COG8  249.23   867.94  113.64
SDHD 1464.84  3508.02    75
SDHD  4.16     3.97       0
SOGA3 81.65   30.76    32.69
SOGA3 0.00    10.77    25.65

These are the Ensembl ids:

SOGA3: ENSG00000214338, ENSG00000255330
SDHD: ENSG00000204370, ENSG00000255292
COG8: ENSG00000213380, ENSG00000272617.

This is not for DE analysis, I want to do subtyping. Should I calculate mean or median on that? How should I go forward with this duplicate genes?

I will use this data for subtyping. And for subtyping there should not be any duplicate genes.

Then you can add them up: if each transcript (ENST..) is of a gene and your subtyping method doesn't handle multiple transcripts, assign them to the same gene. BTW how do you do subtyping (I am looking for applying those methods myself)

just summing them is fine or not should I calculate median or mean on those duplicates?

Then you need to judge if each of these identifiers are really the same gene or not. I would say they are but that's up to you. Why do you think you need to calculate the median or mean? Aren't you interested in knowing how much each gene is being transcripted?

Let's look into this a bit more deeply. For instance:

HUGO: SOGA3
Ensembl 1: ENSG00000214338
Ensembl 2: ENSG00000255330

The Ensembl pages (linked above) for both ENSG00000214338 and ENSG00000255330 show that these genes have 3 paralogs and each lists the other as one of the paralogs. Your other two examples are similar, they are also paralogs of each other. You can also see that the position of each gene in the genome is slightly different:
ENSG00000255330 : Chromosome 6: 127,438,406-127,519,001 reverse strand.
ENSG00000214338 : Chromosome 6: 127,472,794-127,519,191 reverse strand.
So these two genes refer to two distinct regions of the genome. And, as expected, they have different transcripts. However, two of those transcripts (ENST00000481848.6 from ENSG00000255330 and ENST00000525778.5 from ENSG00000214338) code for the same protein sequence. This means that we have one protein but two different genes for it. To illustrate, here are the two transcripts aligned:
ENST00000525778.5       CCGTGTTCGAGCGGTTTTGGGATCCTCTGCCAAGTACGGGATAAGAGAAT
ENST00000481848.6       CTC-----------------------------------------------
                        *

ENST00000525778.5       AGAGCAAGGAGTTTGAGCGCTGCTGCCAATCTTCCATCTCGGGCGTGGCT
ENST00000481848.6       --------------------------------------------------

ENST00000525778.5       CTGGCCTTTTTGTCTCTTTATCCCGCCACTCCCCACCCCCGTCCCCCGTC
ENST00000481848.6       --------------------------------------------------

ENST00000525778.5       CCCCTCCCCGCCTTCTTCTCTTTCCTCTTGAAAGTAACCTCTCTTTGGTT
ENST00000481848.6       -------------------------------------------TTTGGTT
                                                                   *******

ENST00000525778.5       ACTTCTACCCTATCGTCCCTTGCAATCCCGCCTTCTCTCCTGCCCAGTTC
ENST00000481848.6       ACTTCTACCCTATCGTCCCTTGCAATCCCGCCTTCTCTCCTGCCCAGTTC
                        **************************************************

ENST00000525778.5       GGTGAGGTCAGCTCATGAACATCTCTGCACTTCCTGGGGAGCGTCTGGCC
ENST00000481848.6       GGTGAGGTCAGCTCATGAACATCTCTGCACTTCCTGGGGAGCGTCTGGCC
                        **************************************************

ENST00000525778.5       CAGCTTGGTACCTCGCGTGTTAGCCTGCTGAGAGCTGCAGCTCCAGCCCA
ENST00000481848.6       CAGCTTGGTACCTCGCGTGTTAGCCTGCTGAGAGCTGCAGCTCCAGCCCA
                        **************************************************

ENST00000525778.5       GCCGCGGACCGCCGGGAAGTCACCCGCTCTGGGGGCAGCTGAATCTTGAC
ENST00000481848.6       GCCGCGGACCGCCGGGAAGTCACCCGCTCTGGGGGCAGCTGAATCTTGAC
                        **************************************************

ENST00000525778.5       GTTACCGCCATTGTTTGGAGAAATCGAGTGGGTGTGGAACTGCGGCGAGC
ENST00000481848.6       GTTACCGCCATTGTTTGGAGAAATCGAGTGGGTGTGGAACTGCGGCGAGC
                        **************************************************

ENST00000525778.5       GAGCCCAAAGTCATGCTCAGCCGGAAAGTCGCATACAATAATATGAAGTT
ENST00000481848.6       GAGCCCAAAGTCATGCTCAGCCGGAAAGTCGCATACAATAATATGAAGTT
                        **************************************************

ENST00000525778.5       GAACTTTCCTCGGCGCAGAAAGGCACAATAATAAACCCTGAGATCCCCGG
ENST00000481848.6       GAACTTTCCTCGGCGCAGAAAGGCACAATAATAAACCCTGAGATCCCCGG
                        **************************************************

ENST00000525778.5       AATATCAGAAGGAAAGTGCTGTTGATTTGTTCAGCTACTTGGAAACCTCG
ENST00000481848.6       AATATCAGAAGGAAAGTGCTGTTGATTTGTTCAGCTACTTGGAAACCTCG
                        **************************************************

ENST00000525778.5       CCTCTGCTGAGAAGTGAGCAAGTTCTGGATTTCTAAAGGGTTGGGGATGA
ENST00000481848.6       CCTCTGCTGAGAA-------------------------------------
                        *************

ENST00000525778.5       ATGATGGGCTTAGGAAGGAAGGTTGGCGGCGGTGAAGAGGGAAGAGGATG
ENST00000481848.6       --------CTTAGGAAGGAAGGTTGGCGGCGGTGAAGAGGGAAGAGGATG
                                ******************************************

ENST00000525778.5       ATCTAGTAATGGGCTCTTCCAAATCGCATCTGGTAGGTTCATAGCCATGA
ENST00000481848.6       ATCTAGTAATGGGCTCTTCCAAATCGCATCTGGTAGGTTCATAGCCATGA
                        **************************************************

ENST00000525778.5       GTCAGCCACCAATCGGGGGCGCTGCCCCTGCCACAGCAGCAGCTTCTCCC
ENST00000481848.6       GTCAGCCACCAATCGGGGGCGCTGCCCCTGCCACAGCAGCAGCTTCTCCC
                        **************************************************

ENST00000525778.5       GCCGCTGCAGCCACTGAGGCTCGTCTGCACCCGGAGGGCAGCAGCAGAAA
ENST00000481848.6       GCCGCTGCAGCCACTGAGGCTCGTCTGCACCCGGAGGGCAGCAGCAGAAA
                        **************************************************

ENST00000525778.5       GCAGCAGCGAGCTCAGTCACCCGCTAGACCGAGGGACAGTTCGCTGCGAC
ENST00000481848.6       GCAGCAGCGAGCTCAGTCACCCGCTAGACCGAGGGACAGTTCGCTGCGAC
                        **************************************************

ENST00000525778.5       AGACAATCGCGGCCACCCGGTCCCCAGTGGGGGCCGGCACTAAACTCAAT
ENST00000481848.6       AGACAATCGCGGCCACCCGGTCCCCAGTGGGGGCCGGCACTAAACTCAAT
                        **************************************************

ENST00000525778.5       TCTGTTCGGCAGCAGCAGCTGCAGCAGCAGCAACAGCAGGGTAACAAGAC
ENST00000481848.6       TCTGTTCGGCAGCAGCAGCTGCAGCAGCAGCAACAGCAGGGTAACAAGAC
                        **************************************************

ENST00000525778.5       CGGGAGCCGCACGGGGCCGCCCGCCAGCATCCGAGGAGGCGGAGGTGGGG
ENST00000481848.6       CGGGAGCCGCACGGGGCCGCCCGCCAGCATCCGAGGAGGCGGAGGTGGGG
                        **************************************************

ENST00000525778.5       CCGAGAAGGCGACGCCCCTGGCGCCCAAAGGGGCCGCTCCGGGAGCTGTC
ENST00000481848.6       CCGAGAAGGCGACGCCCCTGGCGCCCAAAGGGGCCGCTCCGGGAGCTGTC
                        **************************************************

ENST00000525778.5       CAGCCTGTGGCTGGTGCTGAAGCGGCCCCCGCCGCGACCCTGGCCGCTTT
ENST00000481848.6       CAGCCTGTGGCTGGTGCTGAAGCGGCCCCCGCCGCGACCCTGGCCGCTTT
                        **************************************************

ENST00000525778.5       GGGCGGCAGGAGGCCCGGACCGCCCGAGGAGCCGCCTCGGGAGCTAGAGT
ENST00000481848.6       GGGCGGCAGGAGGCCCGGACCGCCCGAGGAGCCGCCTCGGGAGCTAGAGT
                        **************************************************

ENST00000525778.5       CCGTGCCCTCGAAACTCGGGGAACCCCCTCCGCTCGGAGAAGGGGGAGGA
ENST00000481848.6       CCGTGCCCTCGAAACTCGGGGAACCCCCTCCGCTCGGAGAAGGGGGAGGA
                        **************************************************

ENST00000525778.5       GGGGGCGGGGAAGGAGGAGGAGCCGGCGGCGGCTCCGGGGAAAGGGAGGG
ENST00000481848.6       GGGGGCGGGGAAGGAGGAGGAGCCGGCGGCGGCTCCGGGGAAAGGGAGGG
                        **************************************************

ENST00000525778.5       GGGCGCTCCGCAGCCGCCGCCGCCCAGGGGCTGGCGAGGGAAAGGCGTAC
ENST00000481848.6       GGGCGCTCCGCAGCCGCCGCCGCCCAGGGGCTGGCGAGGGAAAGGCGTAC
                        **************************************************

ENST00000525778.5       GCGCTCAGCAGAGGGGCGGCAGCGGCGGGGAGGGGGCCTCCCCTTCTCCA
ENST00000481848.6       GCGCTCAGCAGAGGGGCGGCAGCGGCGGGGAGGGGGCCTCCCCTTCTCCA
                        **************************************************

ENST00000525778.5       TCCTCCTCTTCTGCGGGCAAAACCCCAGGAACCGGCAGCAGAAACTCCGG
ENST00000481848.6       TCCTCCTCTTCTGCGGGCAAAACCCCAGGAACCGGCAGCAGAAACTCCGG
                        **************************************************

ENST00000525778.5       AAGCGGCGTTGCGGGGGGCGGCAGCGGTGGTGGAGGGAGCTACTGGAAAG
ENST00000481848.6       AAGCGGCGTTGCGGGGGGCGGCAGCGGTGGTGGAGGGAGCTACTGGAAAG
                        **************************************************

ENST00000525778.5       AAGGATGTCTGCAGTCTGAGCTCATCCAGTTCCATCTCAAGAAGGAGCGG
ENST00000481848.6       AAGGATGTCTGCAGTCTGAGCTCATCCAGTTCCATCTCAAGAAGGAGCGG
                        **************************************************

ENST00000525778.5       GCGGCAGCGGCGGCGGCCGCGGCTCAGATGCACGCTAAGAACGGCGGCGG
ENST00000481848.6       GCGGCAGCGGCGGCGGCCGCGGCTCAGATGCACGCTAAGAACGGCGGCGG
                        **************************************************

ENST00000525778.5       CAGCAGTAGCCGCAGCTCCCCGGTGTCTGGCCCCCCTGCCGTTTGCGAGA
ENST00000481848.6       CAGCAGTAGCCGCAGCTCCCCGGTGTCTGGCCCCCCTGCCGTTTGCGAGA
                        **************************************************

ENST00000525778.5       CCCTGGCCGTCGCCTCCGCCTCCCCAATGGCGGCGGCGGCGGAGGGCCCC
ENST00000481848.6       CCCTGGCCGTCGCCTCCGCCTCCCCAATGGCGGCGGCGGCGGAGGGCCCC
                        **************************************************

ENST00000525778.5       CAGCAGAGCGCAGAGGGCAGCGCGAGCGGCGGGGGCATGCAGGCGGCAGC
ENST00000481848.6       CAGCAGAGCGCAGAGGGCAGCGCGAGCGGCGGGGGCATGCAGGCGGCAGC
                        **************************************************

ENST00000525778.5       GCCCCCTTCGTCGCAGCCGCACCCGCAGCAGCTCCAAGAGCAGGAAGAAA
ENST00000481848.6       GCCCCCTTCGTCGCAGCCGCACCCGCAGCAGCTCCAAGAGCAGGAAGAAA
                        **************************************************

ENST00000525778.5       TGCAAGAGGAGATGGAGAAGCTGCGAGAGGAAAACGAGACTCTCAAGAAC
ENST00000481848.6       TGCAAGAGGAGATGGAGAAGCTGCGAGAGGAAAACGAGACTCTCAAGAAC
                        **************************************************

ENST00000525778.5       GAGATCGATGAGCTGAGAACCGAGATGGACGAGATGAGGGACACTTTCTT
ENST00000481848.6       GAGATCGATGAGCTGAGAACCGAGATGGACGAGATGAGGGACACTTTCTT
                        **************************************************

ENST00000525778.5       CGAGGAGGATGCCTGTCAACTGCAGGAAATGCGCCACGAGTTGGAGAGAG
ENST00000481848.6       CGAGGAGGATGCCTGTCAACTGCAGGAAATGCGCCACGAGTTGGAGAGAG
                        **************************************************

ENST00000525778.5       CCAACAAAAACTGCCGGATCCTGCAGTACCGCCTCCGCAAAGCCGAGCGC
ENST00000481848.6       CCAACAAAAACTGCCGGATCCTGCAGTACCGCCTCCGCAAAGCCGAGCGC
                        **************************************************

ENST00000525778.5       AAAAGGCTCCGCTACGCCCAGACCGGGGAAATCGACGGGGAGCTGTTGCG
ENST00000481848.6       AAAAGGCTCCGCTACGCCCAGACCGGGGAAATCGACGGGGAGCTGTTGCG
                        **************************************************

ENST00000525778.5       CAGCCTGGAGCAGGACCTCAAGGTTGCAAAGGATGTATCTGTGAGACTTC
ENST00000481848.6       CAGCCTGGAGCAGGACCTCAAGGTTGCAAAGGATGTATCTGTGAGACTTC
                        **************************************************

ENST00000525778.5       ACCATGAATTAGAAAATGTGGAAGAAAAGAGAACAACAACAGAAGATGAA
ENST00000481848.6       ACCATGAATTAGAAAATGTGGAAGAAAAGAGAACAACAACAGAAGATGAA
                        **************************************************

ENST00000525778.5       AATGAGAAACTGAGGCAACAGCTCATAGAAGTTGAAATTGCAAAGCAAGC
ENST00000481848.6       AATGAGAAACTGAGGCAACAGCTCATAGAAGTTGAAATTGCAAAGCAAGC
                        **************************************************

ENST00000525778.5       TTTACAGAATGAACTGGAAAAAATGAAAGAGTTATCCTTAAAAAGAAGAG
ENST00000481848.6       TTTACAGAATGAACTGGAAAAAATGAAAGAGTTATCCTTAAAAAGAAGAG
                        **************************************************

ENST00000525778.5       GAAGCAAAGATTTGCCAAAATCTGAAAAAAAGGCTCAACAGACTCCCACA
ENST00000481848.6       GAAGCAAAGATTTGCCAAAATCTGAAAAAAAGGCTCAACAGACTCCCACA
                        **************************************************

ENST00000525778.5       GAGGAGGACAATGAAGATCTGAAGTGCCAGCTGCAGTTTGTTAAGGAAGA
ENST00000481848.6       GAGGAGGACAATGAAGATCTGAAGTGCCAGCTGCAGTTTGTTAAGGAAGA
                        **************************************************

ENST00000525778.5       AGCCGCTTTGATGAGAAAGAAAATGGCCAAGATTGATAAAGAAAAGGACA
ENST00000481848.6       AGCCGCTTTGATGAGAAAGAAAATGGCCAAGATTGATAAAGAAAAGGACA
                        **************************************************

ENST00000525778.5       GATTTGAACACGAGCTCCAGAAGTACAGATCCTTTTATGGGGATCTGGAC
ENST00000481848.6       GATTTGAACACGAGCTCCAGAAGTACAGATCCTTTTATGGGGATCTGGAC
                        **************************************************

ENST00000525778.5       AGTCCTTTGCCCAAAGGAGAAGCCGGAGGCCCTCCCAGCACTAGGGAGGC
ENST00000481848.6       AGTCCTTTGCCCAAAGGAGAAGCCGGAGGCCCTCCCAGCACTAGGGAGGC
                        **************************************************

ENST00000525778.5       CGAGCTCAAGCTACGGCTAAGGCTGGTGGAGGAAGAAGCCAACATCCTGG
ENST00000481848.6       CGAGCTCAAGCTACGGCTAAGGCTGGTGGAGGAAGAAGCCAACATCCTGG
                        **************************************************

ENST00000525778.5       GCAGGAAAATCGTCGAACTGGAGGTGGAGAACAGAGGCCTGAAGGCGGAA
ENST00000481848.6       GCAGGAAAATCGTCGAACTGGAGGTGGAGAACAGAGGCCTGAAGGCGGAA
                        **************************************************

ENST00000525778.5       CTGGACGACCTTAGGGGCGATGACTTCAACGGCTCGGCCAACCCGCTCAT
ENST00000481848.6       CTGGACGACCTTAGGGGCGATGACTTCAACGGCTCGGCCAACCCGCTCAT
                        **************************************************

ENST00000525778.5       GAGGGAGCAGAGCGAATCCCTGTCGGAGCTGCGGCAGCACCTGCAGCTGG
ENST00000481848.6       GAGGGAGCAGAGCGAATCCCTGTCGGAGCTGCGGCAGCACCTGCAGCTGG
                        **************************************************

ENST00000525778.5       TGGAAGACGAGACGGAGCTGCTGCGGAGGAACGTGGCCGACCTGGAGGAG
ENST00000481848.6       TGGAAGACGAGACGGAGCTGCTGCGGAGGAACGTGGCCGACCTGGAGGAG
                        **************************************************

ENST00000525778.5       CAGAACAAGCGCATCACGGCGGAGCTCAACAAGTACAAGTACAAGTCCGG
ENST00000481848.6       CAGAACAAGCGCATCACGGCGGAGCTCAACAAGTACAAGTACAAGTCCGG
                        **************************************************

ENST00000525778.5       CGGCCACGACAGCGCGCGGCACCACGACAACGCCAAGACCGAGGCCCTGC
ENST00000481848.6       CGGCCACGACAGCGCGCGGCACCACGACAACGCCAAGACCGAGGCCCTGC
                        **************************************************

ENST00000525778.5       AGGAGGAGCTGAAGGCGGCGCGCCTGCAGATCAACGAGCTCAGCGGCAAG
ENST00000481848.6       AGGAGGAGCTGAAGGCGGCGCGCCTGCAGATCAACGAGCTCAGCGGCAAG
                        **************************************************

ENST00000525778.5       GTCATGCAGCTGCAGTACGAGAACCGCGTGCTTATGTCCAACATGCAGCG
ENST00000481848.6       GTCATGCAGCTGCAGTACGAGAACCGCGTGCTTATGTCCAACATGCAGCG
                        **************************************************

ENST00000525778.5       CTACGACCTGGCCTCGCACCTGGGCATCCGCGGCAGCCCCCGCGACAGCG
ENST00000481848.6       CTACGACCTGGCCTCGCACCTGGGCATCCGCGGCAGCCCCCGCGACAGCG
                        **************************************************

ENST00000525778.5       ACGCCGAGAGCGACGCGGGCAAGAAGGAGAGCGACGACGACTCGCGGCCT
ENST00000481848.6       ACGCCGAGAGCGACGCGGGCAAGAAGGAGAGCGACGACGACTCGCGGCCT
                        **************************************************

ENST00000525778.5       CCGCACCGCAAGCGCGAAGGGCCCATCGGCGGCGAGAGCGACTCGGAGGA
ENST00000481848.6       CCGCACCGCAAGCGCGAAGGGCCCATCGGCGGCGAGAGCGACTCGGAGGA
                        **************************************************

ENST00000525778.5       GGTGCGCAACATCCGCTGCCTCACGCCCACTCGCTCCTTCTACCCGGCGC
ENST00000481848.6       GGTGCGCAACATCCGCTGCCTCACGCCCACTCGCTCCTTCTACCCGGCGC
                        **************************************************

ENST00000525778.5       CCGGGCCCTGGCCCAAGAGCTTCTCCGATCGGCAGCAGATGAAGGACATC
ENST00000481848.6       CCGGGCCCTGGCCCAAGAGCTTCTCCGATCGGCAGCAGATGAAGGACATC
                        **************************************************

ENST00000525778.5       CGCTCGGAGGCCGAGCGCCTGGGCAAGACCATCGACCGGCTCATCGCCGA
ENST00000481848.6       CGCTCGGAGGCCGAGCGCCTGGGCAAGACCATCGACCGGCTCATCGCCGA
                        **************************************************

ENST00000525778.5       CACGAGCACCATCATCACCGAGGCGCGCATCTACGTGGCCAACGGGGACC
ENST00000481848.6       CACGAGCACCATCATCACCGAGGCGCGCATCTACGTGGCCAACGGGGACC
                        **************************************************

ENST00000525778.5       TGTTCGGACTCATGGACGAGGAGGACGACGGCAGCCGCATCCGGGAGCAC
ENST00000481848.6       TGTTCGGACTCATGGACGAGGAGGACGACGGCAGCCGCATCCGGGAGCAC
                        **************************************************

ENST00000525778.5       GAGCTGCTCTACCGCATCAACGCTCAGATGAAGGCCTTCCGCAAGGAGCT
ENST00000481848.6       GAGCTGCTCTACCGCATCAACGCTCAGATGAAGGCCTTCCGCAAGGAGCT
                        **************************************************

ENST00000525778.5       GCAGACCTTCATCGACCGCCTCGAGGTGCCCAAGTCTGCGGACGACCGCG
ENST00000481848.6       GCAGACCTTCATCGACCGCCTCGAGGTGCCCAAGTCTGCGGACGACCGCG
                        **************************************************

ENST00000525778.5       GCGCCGAGGAGCCCATTTCCGTGAGTCAGATGTTCCAGCCTATCATTTTA
ENST00000481848.6       GCGCCGAGGAGCCCATTTCCGTGAGTCAGATGTTCCAGCCTATCATTTTA
                        **************************************************

ENST00000525778.5       CTTATTCTCATTCTTGTATTATTTTCATCACTTTCTTACACAACAATATT
ENST00000481848.6       CTTATTCTCATTCTTGTATTATTTTCATCACTTTCTTACACAACAATATT
                        **************************************************

ENST00000525778.5       TAAACTTGTCTTCCTTTTTACACTGTTTTTTGTACTGTAAATCTTTCATC
ENST00000481848.6       TAAACTTGTCTTCCTTTTTACACTGTTTTTTGTACTGTAAATCTTTCATC
                        **************************************************

ENST00000525778.5       ATTTACCATTCATTGTAGTATTTTCAGTTTGTTTATTTTGTTCACCCTTC
ENST00000481848.6       ATTTACCATTCATTGTAGTATTTTCAGTTTGTTTATTTTGTTCACCCTTC
                        **************************************************

ENST00000525778.5       AAGACAAGAAGTAAAAGAAGTATAATTTCTGTAGTAACCAATGCTATAAA
ENST00000481848.6       AAGACAAGAAGTAAAAGAAGTATAATTTCTGTAGTAACCAATGCTATAAA
                        **************************************************

ENST00000525778.5       AACACTGAAGACTGCTTATTTCTTTAAAAAGATACAACTCATCTTACCAA
ENST00000481848.6       AACACTGAAGACTGCTTATTTCTTTAAAAAGATACAACTCATCTTACCAA
                        **************************************************

ENST00000525778.5       GACCAAA-------------------------------------------
ENST00000481848.6       GACCAAATTCAATAAGAAGCCCAAACACTAAAATATTTCAGGCTTTATTT
                        *******

ENST00000525778.5       --------------------------------------------------
ENST00000481848.6       TAAAGGCAAGTGAGACTGCTTCAAATAAAACAACTTCAAGCTTCCAAGAA

[. . .]

The ENST00000481848.6 transcript continues quite a bit after this but there is no more similarity.
The transcripts are not identical, but are very similar. Now, look at the proteins each transcript encodes:
ENST00000481848.6  MSQPPIGGAAPATAAASPAAAATEARLHPEGSSRKQQRAQSPARPRDSSL
ENST00000525778.5  MSQPPIGGAAPATAAASPAAAATEARLHPEGSSRKQQRAQSPARPRDSSL
                   **************************************************

ENST00000481848.6  RQTIAATRSPVGAGTKLNSVRQQQLQQQQQQGNKTGSRTGPPASIRGGGG
ENST00000525778.5  RQTIAATRSPVGAGTKLNSVRQQQLQQQQQQGNKTGSRTGPPASIRGGGG
                   **************************************************

ENST00000481848.6  GAEKATPLAPKGAAPGAVQPVAGAEAAPAATLAALGGRRPGPPEEPPREL
ENST00000525778.5  GAEKATPLAPKGAAPGAVQPVAGAEAAPAATLAALGGRRPGPPEEPPREL
                   **************************************************

ENST00000481848.6  ESVPSKLGEPPPLGEGGGGGGEGGGAGGGSGEREGGAPQPPPPRGWRGKG
ENST00000525778.5  ESVPSKLGEPPPLGEGGGGGGEGGGAGGGSGEREGGAPQPPPPRGWRGKG
                   **************************************************

ENST00000481848.6  VRAQQRGGSGGEGASPSPSSSSAGKTPGTGSRNSGSGVAGGGSGGGGSYW
ENST00000525778.5  VRAQQRGGSGGEGASPSPSSSSAGKTPGTGSRNSGSGVAGGGSGGGGSYW
                   **************************************************

ENST00000481848.6  KEGCLQSELIQFHLKKERAAAAAAAAQMHAKNGGGSSSRSSPVSGPPAVC
ENST00000525778.5  KEGCLQSELIQFHLKKERAAAAAAAAQMHAKNGGGSSSRSSPVSGPPAVC
                   **************************************************

ENST00000481848.6  ETLAVASASPMAAAAEGPQQSAEGSASGGGMQAAAPPSSQPHPQQLQEQE
ENST00000525778.5  ETLAVASASPMAAAAEGPQQSAEGSASGGGMQAAAPPSSQPHPQQLQEQE
                   **************************************************

ENST00000481848.6  EMQEEMEKLREENETLKNEIDELRTEMDEMRDTFFEEDACQLQEMRHELE
ENST00000525778.5  EMQEEMEKLREENETLKNEIDELRTEMDEMRDTFFEEDACQLQEMRHELE
                   **************************************************

ENST00000481848.6  RANKNCRILQYRLRKAERKRLRYAQTGEIDGELLRSLEQDLKVAKDVSVR
ENST00000525778.5  RANKNCRILQYRLRKAERKRLRYAQTGEIDGELLRSLEQDLKVAKDVSVR
                   **************************************************

ENST00000481848.6  LHHELENVEEKRTTTEDENEKLRQQLIEVEIAKQALQNELEKMKELSLKR
ENST00000525778.5  LHHELENVEEKRTTTEDENEKLRQQLIEVEIAKQALQNELEKMKELSLKR
                   **************************************************

ENST00000481848.6  RGSKDLPKSEKKAQQTPTEEDNEDLKCQLQFVKEEAALMRKKMAKIDKEK
ENST00000525778.5  RGSKDLPKSEKKAQQTPTEEDNEDLKCQLQFVKEEAALMRKKMAKIDKEK
                   **************************************************

ENST00000481848.6  DRFEHELQKYRSFYGDLDSPLPKGEAGGPPSTREAELKLRLRLVEEEANI
ENST00000525778.5  DRFEHELQKYRSFYGDLDSPLPKGEAGGPPSTREAELKLRLRLVEEEANI
                   **************************************************

ENST00000481848.6  LGRKIVELEVENRGLKAELDDLRGDDFNGSANPLMREQSESLSELRQHLQ
ENST00000525778.5  LGRKIVELEVENRGLKAELDDLRGDDFNGSANPLMREQSESLSELRQHLQ
                   **************************************************

ENST00000481848.6  LVEDETELLRRNVADLEEQNKRITAELNKYKYKSGGHDSARHHDNAKTEA
ENST00000525778.5  LVEDETELLRRNVADLEEQNKRITAELNKYKYKSGGHDSARHHDNAKTEA
                   **************************************************

ENST00000481848.6  LQEELKAARLQINELSGKVMQLQYENRVLMSNMQRYDLASHLGIRGSPRD
ENST00000525778.5  LQEELKAARLQINELSGKVMQLQYENRVLMSNMQRYDLASHLGIRGSPRD
                   **************************************************

ENST00000481848.6  SDAESDAGKKESDDDSRPPHRKREGPIGGESDSEEVRNIRCLTPTRSFYP
ENST00000525778.5  SDAESDAGKKESDDDSRPPHRKREGPIGGESDSEEVRNIRCLTPTRSFYP
                   **************************************************

ENST00000481848.6  APGPWPKSFSDRQQMKDIRSEAERLGKTIDRLIADTSTIITEARIYVANG
ENST00000525778.5  APGPWPKSFSDRQQMKDIRSEAERLGKTIDRLIADTSTIITEARIYVANG
                   **************************************************

ENST00000481848.6  DLFGLMDEEDDGSRIREHELLYRINAQMKAFRKELQTFIDRLEVPKSADD
ENST00000525778.5  DLFGLMDEEDDGSRIREHELLYRINAQMKAFRKELQTFIDRLEVPKSADD
                   **************************************************

ENST00000481848.6  RGAEEPISVSQMFQPIILLILILVLFSSLSYTTIFKLVFLFTLFFVL
ENST00000525778.5  RGAEEPISVSQMFQPIILLILILVLFSSLSYTTIFKLVFLFTLFFVL
                   ***********************************************

As you can see, they are 100% identical. So, we have two genes (distinct regions of the DNA) which produce the same protein. Because the gene product is the same, the "function" of the genes is also the same so HUGO treats them as the same gene and gives them the same name.
Whether or not you want to treat them as one gene or two will depend on what you are doing. If you are comparing the relative abundance of each transcript, you must treat them separately. If you want to examine the expression level of the protein product and are using the RNASEM data as a proxy for protein information, then you should treat them as one "gene" and simply sum your RNASEM values.

There is no one-to-one mapping of gene ids from one database to the other. Ensembl (who maintain Ensembl ENSG IDs), ncbi (who maintain EntrezGene IDs and RefSeq transcript ids) and HUGO (who maintain gene symbols/names) all have different ideas about what a gene is, which part of the genome belongs to which gene etc.
In fact, it's worse than that because HUGO gene IDs refer to the abstract idea of a gene, rather than a particular part of any chromosome.
My suggestion would be to leave the IDs as Ensembl IDs (as these actually refer to a physical part of a chromosome, which is what you are quantifying), right up until you have to report the results to someone who really NEEDS gene names. If you still have duplicates, you could sum, or pick the most interesting, but I would suggest just leaving the duplicates in there, and treat it as thought there are two genes with the same name.


text:
SPAdes error during assembly
I try to perform a hybrid assembly using Unicycler. Unicycler use SPAdes to assembly the illumina sequences, and the program crashes while in SPAdes. The spades.log file says something about an invalid argument and error code 22. I have run the spades.py --test, and everything seems to work fine.
Command line: /home/user/SPAdes-3.11.1-Linux/bin/spades.py    -1  /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz -2  /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz -o  /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction --threads   8   --only-error-correction

System information:
  SPAdes version: 3.11.1
  Python version: 2.7.12
  OS: Linux-4.10.0-40-generic-x86_64-with-Ubuntu-16.04-xenial

Output dir: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction
Mode: ONLY read error correction (without assembling)
Debug mode is turned OFF

Dataset parameters:
  Multi-cell mode (you should set '--sc' flag if input data was obtained with MDA (single-cell) technology or --meta flag if processing metagenomic dataset)
  Reads:
    Library number: 1, library type: paired-end
      orientation: fr
      left reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz']
      right reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz']
      interlaced reads: not specified
      single reads: not specified
Read error correction parameters:
  Iterations: 1
  PHRED offset will be auto-detected
  Corrected reads will be compressed (with gzip)
Other parameters:
  Dir for temp files: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/tmp
  Threads: 8
  Memory limit (in Gb): 87

======= SPAdes pipeline started. Log can be found here: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/spades.log

===== Read error correction started.

== Running read error correction tool: /home/user/SPAdes-3.11.1-Linux/bin/hammer /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info

  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  83)   Starting BayesHammer, built from N/A, git revision N/A
  0:00:00.001     4M / 4M    INFO    General                 (main.cpp                  :  84)   Loading config from /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info
  0:00:00.045     4M / 4M    INFO    General                 (memory_limit.hpp          :  51)   Memory limit set to 87 Gb
  0:00:00.045     4M / 4M    INFO    General                 (main.cpp                  :  93)   Trying to determine PHRED offset
  0:00:00.071     4M / 4M    INFO    General                 (main.cpp                  :  99)   Determined value is 33
  0:00:00.071     4M / 4M    INFO    General                 (hammer_tools.cpp          :  36)   Hamming graph threshold tau=1, k=21, subkmer positions = [ 0 10 ]
  0:00:00.072     4M / 4M    INFO    General                 (main.cpp                  : 120)   Size of aux. kmer data 24 bytes
     === ITERATION 0 begins ===
  0:00:00.100     4M / 4M    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 289)   Building kmer index
  0:00:00.101     4M / 4M    INFO    General                 (kmer_index_builder.hpp    : 106)   Splitting kmer instances into 128 buckets using 8 threads. This might take a while.
  0:00:00.101     4M / 4M    INFO    General                 (file_limit.hpp            :  32)   Open file limit set to 1024
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  74)   Memory available for splitting buffers: 3.62484 Gb
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  82)   Using cell size of 524288
  0:00:00.102     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz
  0:00:57.279     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 1338594 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 2152556 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz
  0:09:23.290     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 3442412 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 4305112 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 111)   Total 4305112 reads processed
  0:14:08.726    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 109)   Starting k-mer counting.
  0:20:36.408    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 115)   K-mer counting done. There are 205776690 kmers in total.
  0:20:36.409    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 121)   Merging temporary buckets.
  0:24:43.505    32M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 298)   Building perfect hash indices
  0:25:40.125   100M / 8G    INFO    General                 (kmer_index_builder.hpp    : 137)   Merging final buckets.
  1:11:46.553   100M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 320)   Index built. Total 95426752 bytes occupied (3.70991 bits per kmer).
  1:11:46.598   100M / 8G    INFO   K-mer Counting           (kmer_data.cpp             : 357)   Arranging kmers in hash map order
  1:12:55.597     3G / 8G    INFO    General                 (main.cpp                  : 155)   Clustering Hamming graph.
  1:21:35.596     3G / 8G    INFO    General                 (main.cpp                  : 162)   Extracting clusters
=== Stack Trace ===
[0x407a6a]
[0x409561]
[0x40c79b]
[0x40129a]
[0x5442b0]
[0x40620d]
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
hammer: /spades/src/common/io/kmers/mmapped_writer.hpp:94: void MMappedWriter::reserve(size_t): Assertion `(intptr_t) MappedRegion != -1L' failed.

== Error ==  system call for: "['/home/user/SPAdes-3.11.1-Linux/bin/hammer', '/home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info']" finished abnormally, err code: -6

In case you have troubles running SPAdes, you can write to spades.support@cab.spbu.ru
or report an issue on our GitHub repository github.com/ablab/spades
Please provide us with params.txt and spades.log files from the output directory.

Answer from @mjoppich, converted from comment:
Late answer, but it looks like you are running out of memory, hence the intptr is invalid.


text:
Proteins that cannot form biofilm?
I am trying to build a machine learning training set for bacterial protein sequences that form biofilm, and that cannot. I collected the positive sequences from the GO ontology website but for negative sequences I am not sure which sequences to incorporate into my training set.
Can anyone point me to resources for proteins sequences that are known to be not capable of forming biofilms?

Maybe I am wrong but all the proteins you don't have as your positive set is your negative set. How do you select the biofilm capable proteins? By those annotated with a specific GO term?

One option is to treat all proteins that are not explicitly known to form a biofilm as proteins that are incapable of doing so, but this will likely result in many false negatives.
You could take a look at one-class classification and positive-unlabeled (PU) learning. Both are techniques specifically designed for these kind of problems where only one of the two classes is known. (In this case: only biofilm forming proteins are known.)


text:
Comparison of gene set enrichment statistics
I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.
As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.
How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply.
I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:
library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")

# Create a design matrix and
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)),
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
##
## $correlation
## [1] -0.02036983
##
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
##
## $correlation
## [1] 0.3024643
##
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
##
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance.
As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:
               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value?

I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications.

Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used inter.gene.cor = NA in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower.
Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (cameraPR) because it doesn't have the same results as in camera and it doesn't allow to estimate the intercorrelation.

Noticed an error in your code - in the call to topTable, you need to include 'sort.by = "none"' or else it reorders by significance, making the 1:40 and 41:80 indices for geneSetTest no longer applicable.The geneSetTest p-values for the three gene sets now come out to 1, 9.04e-13, and 0.97. It seems that the geneSetTest p-value is arguably more "correct" here, since only GS2 truly shows differential expression between samples 1-5 and 6-10. CAMERA says they're all significant!

@NuclearWang Nice catch, that's why I usually use ids2indeces! I am not sure if this is the real case, I'll try to expand the answer using the ALL dataset in a couple of days. Hope this helps! If you find a better answer, please let me know, I am interested in those methods :)

Why is CAMERA giving a significant p-value for the first comparison?

@flies There is only one comparison. Do you refer to the enrichment in the GS2?

@Llopis I refer to GS1. AFAICT, rows 1 through 40 are just a uniform normal distribution across samples, so the test should be negative.

@flies Yes, you are right. The GSEABenchmarkeR aims to find a consensus. My understanding is that the low correlation and the VIF makes strange things. But it might be an effect of a small dataset...

@Llopis thanks for your response on this old question! It's really odd and rather disturbing that this is what you found. I will say that I recently had a case where I used cameraPR where the ids2indices was scrambled with respect to the expression input (I used t-statistics output by limma), and the significance was all negative. Once I fixed the bug that scrambled the indices, I got good results that made biological sense. So camera will give negative results in some cases - at least with fixed inter.gene.cor!


text:
Coordinates to Genome R Studio Script?
I have a very, very long list of genome coordinates of areas. I need to gather a sequence in order to compare it to another database of sequences that I have. How do I produce the genome sequences from these coordinates in R?
I am very new to this and need a detailed answer.

What genome? Namely, is it one for which there's already a BSgenome package?

The genome was given to me, it is not on there. My problem, more specifically I suppose, is that I was given a couple files. The genome was run through RepeatMasker to obtain coordinates of transposable elements. So, I have the coordinates from the map file that was produced from RepeatMasker and I would like to identify the sequences of these coordinates from my genome (which is too big to view). After I obtain the sequences, I will need to compare those sequences to a library of known elements that was also given to me in another file. So, I am lost especially since I am very new to this

OK, does this really have to be done in R? This particular task is a bit better done on the command line than in R.

Do you know what type of commands this would consist of?

This doesn't use R, but a simple and quite convenient method is to use bedtools getfasta with your fasta file. You'll need to make a BED file out of the RepeatMasker output:
awk 'BEGIN{OFS="\t"}{name=sprintf("%s_%s", $10, $11); strand="+"; if($9 == "C"){strand="-"}; print $5,$6-1,$7,name,".", strand}' repeatmasker.output > rmsk.bed

Then, given that and your genome.fa fasta:
bedtools getfasta -name -s -fi genome.fa -bed rmsk.bed -fo output.fa

This will produce a multiline fasta file names like AluSx_SINE/Alu (due to the -name option) where the sequences are reverse complemented if appropriate.
You can install bedtools either via conda (conda install -c bioconda bedtools) or github.


text:
Can we use non base-called fast5 files in poretools?
I run the MinION MinKNOW without the live base-calling option. We know there is Metrichor and Albacore to perform base-calling after this process. However, I have not done any base-calling yet.
My question is: Is it possible to use directly the fast5 files with poretools to extract a fastq file without doing any previous base-calling process? I tried it and I get an empty fastq file.
My reads directory only contains the fast5 folder and I run:
poretools fastq fast5/ > output.fastq
Any ideas why do I get an empty file? What is the difference between base-called fast5 and non-base-called fast5 files?
Now I am trying to do base-calling with Albacore to see if I get a fastq file.

You’ve edited the title - I should add that poretools can work with non base-called FAST5 files with some of its options, but not the fastq option. Some of the metadata options will function.

No, poretools does not do basecalling. The poretools fastq command can be used to extract the FASTQ information from the basecalled FAST5 file (via MinKNOW live-basecalling or albacore). Alternatively, both of these basecallers can export a FASTQ file directly if desired.
The difference between a basecalled FAST5 and a nonbasecalled FAST5 is the presence of the FASTQ within the file. This can be directly viewed using an HDF5 viewer, like h5dump from the commandline, or HDFView from the GUI.

Albacore can be used to create a fastq, directly from your fast5 files. As such there is no more reason to use poretools, expect perhaps for getting some metadata. With regards to metadata I would like to add that the sequencing_summary.txt generated by albacore also contains a lot of metadata.