text:
Plot to show the expression of genes between tumor and normal
I have RNA-seq raw counts data for 50 samples. 20 Normal and 30 tumor. After differential analysis I got 30 DEGs. I want to make a violin plot showing the expression of each gene. I transformed counts to logCPM.
counts:
Genes   Tumor1  Tumor2  Normal1 Normal2
RP11-351J23.1   0   5   6   0
MIR17HG 989 896 8   0
RP11-563N12.2   0   0   0   0
LINC04393   0   2   16  0
RP11-336A10.4   0   0   0   0
DRH6-AB1    53  13  39  9
RP11-115J16.1   0   0   50  6
LINC70518   2   65  0   0

logCPM <- cpm(counts, prior.count=2, log = TRUE)

                Tumor1   Tumor2  Normal1  Normal2
RP11-351J23.1 11.84477 13.09301 15.71337 11.84477
MIR17HG      19.84847 19.79600 16.10350 11.84477
RP11-563N12.2 11.84477 11.84477 11.84477 11.84477
LINC04393     11.84477 12.47723 17.06532 11.84477
RP11-336A10.4 11.84477 11.84477 11.84477 11.84477
DRH6-AB1      15.72257 14.03896 18.32772 19.19283
RP11-115J16.1 11.84477 11.84477 18.68262 18.61229
LINC70518     12.44599 16.08368 11.84477 11.84477

I want to make violin plot for each gene showing expression for Tumors and Normals with this.
I first transformed logCPM. then tried like below:
ggplot2.violinplot(data=logCPM, xName='RP11-351J23.1',
    groupName='Tumor', position=position_dodge(0.8),
    backgroundColor="white", groupColors=c('#999999','#E69F00'),
    legendPosition="top")

From what package is this ggplot2.violinplot function? Could you upload the image you get and why it doesn't meet your requirements?

"easyGgplot2" is the package. I didn't get any plot. something went wrong and no error also.

Could you post your sessionInfo then? And is there any reason why you use this function and not the base ggplot2 functions?

You'll want to use long-form data for everything:
library(dplyr)
library(ggplot2)
logCPM$gene = row.names(logCPM)
d = logCPM %>% gather(Sample, logCPM, -gene)
d$group = c(rep("Tumor", 2), rep("Normal", 2))
geneOfInterest = d %>% filter(gene == 'RP11-351J23.1')
ggplot(geneOfInterest, aes(x=group, y=logCPM)) + geom_violin()

That's a simple example, you can tweak it to meet your needs.

With logCPM mentioned in my question I have used your code. But there're some warnings and error. logCPM$gene = row.names(logCPM) Warning message: In logCPM$gene = row.names(logCPM) : Coercing LHS to a list; d = logCPM %>% gather(Sample, logCPM, -gene) Error in UseMethod("gather_") :    no applicable method for 'gather_' applied to an object of class "list"

Make sure logCPM is a data.frame

yes. but In the plot I can see only Tumor. There is no Normal. https://imgur.com/a/LnxG4DN And how to add asterisks to denote a significant difference?

You somehow filtered out the normal samples. You can figure out how that happened. You can add characters with geom_text(), please see the ggplot documentation.

I followed the code. But dataframe "d" looks like this;                                                gene  Sample   logCPM  group
1  RP11-351J23.1  Tumor1 11.84477  Tumor
2       MIR137HG  Tumor1 19.84847  Tumor
3  RP11-563N12.2  Tumor1 11.84477 Normal
4      LINC00393  Tumor1 11.84477 Normal
5  RP11-336A10.4  Tumor1 11.84477  Tumor
6       DLX6-AS1  Tumor1 15.72257  Tumor
7  RP11-115J16.1  Tumor1 11.84477 Normal
8      LINC00518  Tumor1 12.44599 Normal

Ah, you'll need to change the d$gene assignment line.


text:
Detect gene duplication events from a whole genome in FASTA format
Are there any tools for finding duplications from a whole genome in FASTA format?
Let's say the genome has gene A duplicated somewhere in the genome and gene B triplicated in tandem elsewhere in the genome. Is there a tool that will output a list of gene A and gene B with location and number of copies?
I am currently using phobos to annotate tandem repeats in the genome of C. elegans but it becomes very slow beyond repeats of 100 bp. Does anybody know of a tool for gene-size/larger repeats?
I have tried to use TRAL but I have got some installation and usage issues, described here.
I cannot use the annotation because I would like to apply this tool to a non-reference version of a C. elegans genome, for which no annotation is available.
Specifically, I would like to annotate the Hawaiian strain, which is a strain very divergent to the reference N2 strain. I am not looking into particular genes, just natural genetic variation.
I have not tried BLAST. It would be a possibility but I wonder if there are others.
Cheers.

If you're working on C. elegans isn't this information already available? Surely the duplicated genes will be well known in an organism this well studied!

I would like to do this analysis on a non-reference version of C. elegans.

Please [edit] your question and explain more clearly. If this is for a C. elegans genome, why would you expect it to be significantly different? Can you at least assume it will be similar enough that you will have a known set of genes to look for? Do you only want to look for gene duplications or large repeats in general? Have you tried simply blasting (tBLASTn) all known C. elegans proteins against the genome?

@terdon I hope my question is clearer now. I have not tried BLAST. I am considering it.

@terdon The orthologs will be annotated for C. elegans (e.g. OMA, OrthoDB, etc), but the number or tandem repeats likely will not be annotated.

@charlesdarwin I'm not aware of any tools that output both number of repeats and homology information. That would be a useful DB. I'm hoping someone answers this in the affirmative!

source: Stack Exchange

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

sort: votes


text:
How to merge transcript sequence with same name in a FASTA file.
Suppose we have a fasta file like
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCAT
>Seq1
AAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

I want to produce
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCATAAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

Certainly this python script is not the most elegant solution but it will give you the desired result:
import sys

seqs = {}
with open(sys.argv[1], "r") as fh:
    curr = ""
    for line in fh:
        if line.startswith(">"):
            curr = line
            if curr not in seqs:
                seqs[curr] = ""
        else:
            seqs[curr] += line.strip()

for ident, seq in seqs.items():
    print("{}{}".format(ident, seq))

If you saved the script in seq_merge.py just run it with:
python seq_merge.py your_file.fasta > result.fasta

This is a scenario where awk excels and yields a much simpler solution than the other examples here (though the logic is essentially equivalent to terdon’s Perl solution):
awk '
    />/  { id = $0 }
    !/>/ { seq[id] = seq[id] $0 }
    END  { for (id in seq) print id "\n" seq[id] }
' input.fasta \
> output.fasta

That said, for anything beyond a quick’n’dirty application I strongly recommend using a proper FASTA parser, not an ad-hoc solution.

assuming two lines per fasta record (or you can linearize) and using datamash:
$ cat jeter.fa | paste - - | sort -t $'\t' -k1,1 | datamash -g 1 collapse 2 | tr -d ',' | tr "\t" "\n"

>Seq1
AAAAACGAAACTGTATCCCGTGTTTGTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCAT
>Seq2
CGTGTTTAGCAAAGAAAT

If you don't care about the order of the sequences and your file is small enough to fit in RAM, you can do:
$ perl -nle 'if(/^>/){
                $n=$_
             }
             else{
                $seq{$n}.=$_
            }
            END{
                print "$_\n$seq{$_}" for keys(%seq)
            }' file.fa
>Seq1
GTTGAGAGGTGTATGGACACGAAAAACGAAACTGTATCCCGTGTTTAGCAAAGAAATCATAAAAACGAAACTGTATCCCGTGTTT
>Seq2
CGTGTTTAGCAAAGAAAT

Or, if you are into the whole brevity thing1 :
perl -nle '/^>/?($n=$_):($s{$n}.=$_);}{print"$_\n$s{$_}"for keys%s' file.fa

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


text:
Plotting Ramachandran Plot from more than one PDB file
I want to plot a Ramachandran Plot for around 5000 PDB files in python that is saved in a folder in my system. I am searching for a way to make it without accessing each file by hard-coding each file name. Please tell if there is a way around this.

Is your problem only how to avoid hardcoding file names? If so, please show us your code so we can adapt it. Or do you also need help in making the plot? What do you have so far?

@terdon It seems the flag for homework is missing. No homework policy yet? This is pretty standard introductory task.

@Mithoron no, we have no homework policy. And I hope we don't get one, to be honest. Questions can be good and clear, or bad and vague irrespective of whether they're homework. But feel very free to bring this up in meta if you think we should have one!

You could probably use PyRAMA for that. The only problem is that by PyRAMA plots to display, while you probably would like to save to a file.
One way to get around this is to:

Patch PyRAMA to allow ramachandran function to accept output image filename.
def ramachandran(file_name_list, output_image=None):
    ...
    plt.tight_layout()
    if output_image is not None:
        plt.savefig(output_file, dpi=300)
    else:
        plt.show()

If this works as expected I would really encourage you to submit a pull-request (or at least open an issue) on github.
Now you can write your main script:
from glob import iglob

from rama_lib.ramachandran_calc import ramachandran

for pdb_filename in iglob('/my/dir/with/pdb/', '*.pdb'):
    img_filename = pdb_filename.rsplit('.', 1)[0] + '.png'
    ramachandran((pdb_filename,), img_filename)

I did not test the code, but I hope you get an idea.
There are also other ways to make a ramachandran plot in python. Please check: 1, 2, 3.


text:
biological fluorescence supervised machine learing
Some general questions...
So it appears with most classification learners one must come up with a series of quantifiable variables to associate with each observation.  This might include mean and std deviation intensity values, texture, and zernike polynomials.  Does anyone have any other suggestions as to how one might quantify 3-dimensional fluorescence pixel data features?
Once appropriate quantification's have been acquired for a given class one can then train a classifier.  My next question is - does one then have to segment test images and create a matrix of the same feature variables from that segmented image to use the classifier? Or can one use the raw image itself?
What supervised machine learning techniques will evaluate the raw image data itself - either for training or testing? Can manually annotated data sets or clustered data sets provide a direct model for training data?  That is to say, are there supervised machine learning algorithms which accept pure pixel data as models for training?
Are there reliable pre-trained classifiers that are already implemented for these type of fluorescence data sets?
Any resources or advice is greatly appreciated.
Thanks in advance for any input!

Whether to quantify things at all with depend on the technique you want to use. An end-to-end CNN (convolutional neural network) for example would just be fed the 3D image for training and predictions. However the major downside to CNNs is that you need a LOT of data to train them, though perhaps you can get lucky and find an existing CNN that you can use and just retrain the last couple layers.
Other than that, you're going to need to come up with as many (hopefully useful) metrics as you can. There are then a lot of techniques one could try, ranging from simple things like regression to random forests to SVMs, among many others. In general I think simple regression and random forests are the easiest to get started with.

Excellent.  Thank you!  Now If i understand correctly, one feeds the image in 1 dimension essentially?  That is to say the total number of elements reshaped into one dimension?   From what I've seen the layers are all 1 dimensional unless the images have color layers

Typically yes, though in reality it's a contiguous stretch of memory already, so you're just ignoring the dimensional aspects when running the CNN.

Adding to Devon's correct answer: Depending on the amount of images you got, there is also the possibility to avoid the feature engineering by building a two level architecture:

The first level identifies the regions of interest in your image and scales these regions to the same size.
In the second level, the recognized areas are fed to a (multilayer) neural network/deep neural network to be classified.

By using this approach, the number of needed datapoints is lower compared to a plain CNN method.
The drawback is that identifying the regions of interest itself might be a difficult task as you have either to label these regions manually to apply a machine learning approach at this first level, or you have to "engineer" criteria which make the extraction possible.
I have no idea how these images look like, but maybe you are able to extract the regions using OpenCV.

Very interesting.  I will check this out also.  I have been working with some clustering to isolate some examples of objects of interest.   I haven't perfected this process yet as the fluorescence data is somewhat noisy.
The same question for you as Devon, is it possible to send the image in as a 3d stack array?  Or will the image necessarily be sent in on layer of the stack at a time?  My interpretation thus far is the latter

At least to OpenCV and a CNN you can present the images one at a time in a 3D format, i.e.: width x height x color. That means you total dataset has to have the format number_of_images x width x height x color

sounds good thank you the assistance and advice!

What supervised machine learning techniques will evaluate the raw
  image data itself - either for training or testing?

You can use convolutional neural networks (CNN). They will directly use the raw image data itself.
After applying CNNs, you can also use recurrent neural layers (RNN), to achieve better results (e.g. gated recurrent layers).

Hope this helps. As a reference you have this article:

Gudla PR1,2, Nakayama K2,3, Pegoraro G1,2, Misteli T2. Article:
  SpotLearn: Convolutional Neural Network for Detection of Fluorescence
  In Situ Hybridization (FISH) Signals in High-Throughput Imaging
  Approaches. Send to Cold Spring Harb Symp Quant Biol. 2017;82:57-70.
  doi: 10.1101/sqb.2017.82.033761. Epub 2017 Nov 28.

Please include some links to resources on CNN and RNN

Gudla PR1,2, Nakayama K2,3, Pegoraro G1,2, Misteli T2.
SpotLearn: Convolutional Neural Network for Detection of Fluorescence In Situ Hybridization (FISH) Signals in High-Throughput Imaging Approaches.
Send to
Cold Spring Harb Symp Quant Biol. 2017;82:57-70. doi: 10.1101/sqb.2017.82.033761. Epub 2017 Nov 28.
Abstract and full text link can be found here:
https://www.ncbi.nlm.nih.gov/pubmed/?term=3D+CNN+fluorescence

Please edit your answer to include those links

Authors: Gudla PR1,2, Nakayama K2,3, Pegoraro G1,2, Misteli T2.
Article: SpotLearn: Convolutional Neural Network for Detection of Fluorescence In Situ Hybridization (FISH) Signals in High-Throughput Imaging Approaches. Send to Cold Spring Harb Symp Quant Biol. 2017;82:57-70.
doi: 10.1101/sqb.2017.82.033761. Epub 2017 Nov 28.
you can find this article in https://www.ncbi.nlm.nih.gov/pubmed/
A detailled link is here:
Abstract and full text link can be found here: ncbi.nlm.nih.gov/pubmed/?term=3D+CNN+fluorescence


text:
Comparing two files with overlapping regions and get their associated information
I have two file that have overlapping regions. One file has chromosomes, position, CNV data and samples.Other file has chromosomes, positions and ensemble ids of genes.
Here is some lines of two files.
contig  start   end     type  bi-allelic  S1  S2  S3
chr1    12048   177033    DUP FALSE   -1  22  24
chr1    12048   89237     DUP FALSE   12  10  -1
chr1    2985712 3355300   DEL   TRUE  2   2   2

Second file is like this
chr1    901877  911245  gene_id "ENSG00000187583.6"
chr1    22013   41670 gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

I want output like this.
chr1    22013   41670 gene_id "ENSG00000160075.10"  DUP FALSE   12  10  -1
chr1    2985732 3355185 gene_id "ENSG00000142611.12" DEL   TRUE  2   2   2

Is there any way to do this?

You can use bedtools intersect and cut the results:
$ cat file1.bed
chr1    12048   177033  DUP FALSE   -1  22  24
chr1    12048   89237   DUP FALSE   12  10  -1
chr1    2985712 3355300 DEL TRUE    2   2   2

$ cat file2.gtf
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    22013   41670   gene_id "ENSG00000160075.10"
chr1    2985732 3355185 gene_id "ENSG00000142611.12"

$ bedtools intersect -a file2.gtf -b file1.bed -wb | cut -f 1-4,8-11
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   12  10
chr1    22013   41670   gene_id "ENSG00000160075.10"    DUP FALSE   -1  22
chr1    2985732 3355185 gene_id "ENSG00000142611.12"    DEL TRUE    2   2


text:
How to visualize genome track of gene in specific cell-lines?
I'm trying to make a plot showing genome tracks of specific genes in specific cell-lines of RNA-seq and Chip-seq data.
It should look something like this I have recently seen this encode, but in this I didn't find any MDAMB231 cell-lines.
May I know how to make such plots.
I got the Counts data from CCLE portal.
             BT-549 MDAMB231
RP11-554I8.2    0   246
LINC00707      50   107
   RBM17    17036   10819
  PFKFB3    18098   7668
RP11-393K12.2   0   0

You can use screen.encodeproject.org
It collected a lot of data from multiple cell lines or tissues. And you can export tracks to UCSC browser or other browsers.

Answer from @devon-ryan, converted from comments:
You'll need BAM or bigWig files for each of the tracks. If you can't find any you'll have to make them yourself from public fastq files. The fastq files can be obtained from ENA, SRA, or any of the other standard databases.

Like others have mentioned, these plots come from BAM, BED, and BigWig files.
Those files contain the data that you can give to tools like the UCSC Genome Browser, IGV,  HiGlass/Gosling, or other bioinformatics visualization software.
For the specific cell line you mention, you'll have to find the relevant RNA-seq and ChIP-seq data from some data repository, like CCLE, ENA, SRA, or GEO.
You may be to put in some work converting available files to one of these compatible formats before plotting them.


text:
Annotation with Prokka or RAST
I was experimenting Prokka and RAST annotation tools. So, I took a well-annotated swinepox virus genome from genebank (NCBI Reference Sequence: NC_003389.1).
I ran those sequences on Prokka and RAST Seed server at the same time. I can see that only a few (may be around 1%) of the genes were annotated. Most of them were predicted as hypothetical protein. And the results were comparable between Prokka and RAST.
I would assume that these tools look for similar sequences in NCBI and find the best-match protein. But looks like that is not the case.
I am looking for prediction tool(s) that will:

find a well annotated swinepox virus genome in Genbank to predict most of the proteins.
improve on the prediction of Prokka and Rast to assign functions to otherwise hypothetical protein

Would prediction tools such as Genemark or other genome annotation tools be beneficial.
Alternatively have I misunderstood the concept of annotation?
Assistance would be very welcome.
I have attached the image for comparison. Left one the Swinepox genome in .gb format and right one is the same genome annotated with Prokka.

Please, provide links to the annotation servers along with all options you've specified.

I used Galaxy server.  I chose viruses under kingdom. All the other parameters were kept as default.

Prokka works well on bacterial genomes. I quesiton whether it will be useful here.

Both of the annotation tools you used are designed for prokaryotic genome annotation. I would not expect them to work very well for viruses.

You can choose between Bacteria, Archea or Viruses when using prokka. So, I think it should work for viruses too.

@LRJoshi viruses (especially small ones) tend to have overlapping and nested ORFs: something a standard prokaryotic/eukaryotic gene detector (e.g. Glimmer, Genemark, Prodigal) is unlikely to account for.

Prokka can do viruses, but you may want to change which translation table it uses to get better blast results

I was able to get good annotation using prokka by using following commands in prokka.
prokka --proteins reference.gb --outdir annotation  --prefix myprotein contigs.fa

All the reading frames were annotated. Actually the ones that are hypothetical in reference genomes are also annotated.

I definitely agree that Prokka could be repurposed and think this is a good response.

If you are looking for another tool, you could try eggnog mapper: it is an annotation tool based on orthology assignments. The EggNOG database happens to have a section for viral Orthologous Groups so you might obtain better results with it.
Note that emapper annotates coding sequences, not whole genomes. You need an extra step to predict CDS, for instance using Prodigal.
Basic usage with the optimized database for viral models:
prodigal -i my.genome.fna -o my.genes -a my.proteins.faa

then
python emapper.py -i my.proteins.faa --output polb_viruses  -d viruses

Interesting, do you have a source on this? The documentation of eggnog-mapper mentions "the annotation of novel genomes". I am currently using it to annotate MAGs, so I'd rather know if I am doing something wrong.

Eggnog resulted in empty file. Maybe it doesn't work for individual genomes as mentioned by @EliKorvigo

eggnog-mapper takes predicted proteins as input, not whole genomes. You should not use it on raw metagenomes either. You need an extra step to predict CDS with the tool of your choice (e.g. Prodigal). prodigal -i my.genome.fna -o my.genes -a my.proteins.faa then feed my.proteins.faa to eggnog-mapper.

My wording has been quite a bit too harsh, I guess. I have only applied eggnog to metagenomes and metatranscriptomes, because in these cases it's not practical to use group-specific homology databases and HMMs. Eggnog's authors have even emphasised this point in their paper. On the other hand, when you do have a clean assembly, you can use a hierarchical group-specific annotator, such as Prokka. Of course, it doesn't mean one can't use eggnog for genomes, hence I find it necessary to apologise for my previous statement.


text:
Variants in Exome that affect miRNAs related to Schizophrenia
I have 3 exomes results (VCF files), I would like to study variants (SNP, INDEL) that can affect the miRNAs, which are related to a complex disease (Schizophrenia). What is the best way to assess the effects of variants on miRNAs and how can I link between the damaged miRNAs and Schizophrenia. Can I find a help or advice please?

I assume that you want a bioinformatician approach, not an experimental one, and that depends on the data you have (or you can obtain). For instance, how do you measure the damaged miRNAs? What data do you already have?

Do you want effects on miRNAs, or effects on their targets?

I would run them through the Ensembl VEP, which will tell you which genes they hit and the effect they have on them. Then I would filter by consequence "mature_miRNA_variant", which will find you all the variants in your file that fall in the mature miRNAs. Or you could filter instead by the gene biotype "miRNA" which will tell you if the variants fall anywhere in the miRNAs (which will include the full stem loop structure). I would then export the list of miRNAs and try to find out about their targets using miRbase.


text:
ATACseq pipeline report generation issue
I am trying to run ATACseq pipeline:
https://www.bioconductor.org/packages/release/bioc/html/esATAC.html
On single end ATACseq experiments.
conclusion_micro12 <- atacPipe("mm10", fastqInput1 = paste0(base_folder, "Micro_atac12.fastq.gz"),
                   adapter1 = "CAAGCAGAAGACGGCATACGAGATAGCGTAGCGTCTCGTGGGCTCGGAGATGT",
                   threads = 4,
                   createReport = TRUE)

It generated several folders during its run: esATAC_report, esATAC_result, intermediate_results, refdir. However, only the last two gets populated, the others being empty. I am not sure what is going on. I looked over various log files in the intermediate_results folder, but it does not say much. Where should I look to figure out what went wrong? I would expect a report to be generated if everything is fine since I specified createReport = TRUE. While running the pipeline it does not show any errors in Rstudio console.


text:
How to scale the size of heat map and row names font size?
I have an expression data matrix (120X15; 15 samples and 120 genes), my heatmap looks blurred and raw names (gene names) looks very small and can not read. How can I improve my scripts?
Here is the example data
df<-structure(list(X_T0 = c(1.32606106633333, 0.222768019,  1.27419565966667, 0.293489379),
X_T1 = c(1.499956424, 0.293117217, 1.11955345266667, 0.324453864666667),
X_T2 = c(1.118528738, 0.548228048, 0.936319105, 0.846378238666667),
X_T3 = c(1.025082136, 0.841192647333333, 0.957743318666667, 1.10810667833333),
Y_T0 = c(0.424537206666667, 3.144197864, 1.01435004766667, 2.195831453),
Y_T1 = c(0.723243112666667, 0.576764958333333, 0.604810037, 0.896297411666667),
Y_T2 = c(0.335509156333333, 1.32037215366667, 1.209961935, 0.290895345333333),
Y_T3 = c(0.328275209, 1.15039119233333, 1.409758678, 0.0922298696666667),
Z_T0 = c(0.788329993666667, 1.03539976366667, 1.10523299, 0.885098106666667),
Z_T1 = c(1.125292329, 1.00032109266667, 0.877359233, 1.385935293),
Z_T2 = c(2.357924224, 0.740699933666667, 0.585095342, 1.98251017733333),
Z_T3 = c(0.678921448, 0.687992671666667, 1.005347728, 1.05829030233333)), row.names = c("AA", "PQ", "UBQ", "XY"), class = "data.frame")

Script being used
col <- colorRampPalette(c("red","yellow","darkgreen"))(30)
png("S7.png")
heatmap.2(as.matrix(df), Rowv = T, Colv = FALSE, dendrogram = "row", #scale = "row", col = col, density.info = "none", trace = "none", margins = c(7, 15) )
dev.off()

For sample data, this script generates perfect heatmap, but when I have more raw data(genes), this problem comes
I've tried to change these margins = c(7, 15) value, doesn't work well and not sure how to proceed.
Also, I want to reduce the height of the color key

I'd appreciate your help with this.

heatmap.2 is very configurable, and has options to adjust the things you want to fix:
cexRow: changes the size of the row label font.
keysize: numeric value indicating the size of the key.
The size of the key is also affected by the layout of the plot. heatmap.2 splits your plotting device into 4 panes (see the picture below), and you can control the size of the key partly by controlling the size and layout of these 4 panes.

The layout and scale of these panes can be controlled using the lmat, lwid and lhei arguments. For example, in your case, because you're not displaying the column dendrogram in your plot, you could set the height of the first row of the layout to be smaller, which would reduce the size of the colour key: lhei=c(2, 10) (this is just an example, you'd need to experiment to find values which worked well for your specific heatmap).

Make the image itself bigger (e.g., png("S7.png", width=1000, height=1000)).
Having said that, rethink the utility of having the labels there. Are you really going to go through and look at them if you have hundreds of lines? Probably not.


text:
How to interpret PCA output statistically and biologically?
How can I interpret the PCA results statistically for biological data?
I have used FactoMineR and factoextra libraries for PCA
Scripts used:
library(FactoMineR)
res.PCA = PCA(df, scale.unit=TRUE, ncp=4, graph=F )
par(mfrow=c(1,2))
plot.PCA(res.PCA, axes=c(1, 2), choix="ind")
plot.PCA(res.PCA, axes=c(1, 2), choix="var")
dimdesc(res.PCA, axes=c(1,2))

library("factoextra")
fviz_pca_var(res.PCA, arrowsize = 1, labelsize = 3, repel = TRUE, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Here is my PCA output of differnt KO genotypes); how can I explain/interpret the PCA output statistically and if possible biologically: esp, how this colors /contrib explain

Oncidium, did you check it with the wikipedia page on pca or with the help page of the function? Also what kind of data have you used as input (this would help with the biological intepretation)?

@Llopis, This PCA is from targeted proteomics data, sorry I have tagges rna-seq, it should be proteomics.

All that plots like this are telling you is that there are some genes that contribute more to the variability seen between your various samples than others. In an ideal world these genes will also be differentially expressed between your groups, but since we don't live in an ideal world that might not be the case. In general, the longer the line and the closer it is to the X axis the more it contributes to the overall variation (since it's contributing more the PC1, which is 53.5% of the variation in your data).
Since you've asked in the comments for a brief "results and description" section style summary of such a PCA I'll include something brief below. Note that one should not over-interpret PCA plots. The axes are (mathematically) orthogonal and don't necessarily correspond to anything biologically relevant.

We performed a PCA on the variance-stabilized counts to check for batch effects and overall clustering of the data. As can be seen, the "3T" and "5T" groups cluster together along the first principal component, while the "0T" and "1T" samples cluster on the opposite side. This suggests that these groups are more similar to each other, as further evidenced by ...

Assuming you had used the vst() function in DESeq2 (I note you've tagged your post with RNAseq, so perhaps you're using DESeq2) then something like that above would be a reasonable description.

Thanks Dr. Ryan, May I ask you if I can have a example of clear description for PCA results.  I mean how can I put this kind of results in to words when writing an R&D. I have never seen clear interpretation of PCA results in papers I have gone through.

This ends up being rather specific to the topic. In general PCAs are just used for QC to see if there's something odd going on (e.g., an outlier sample or things clustering oddly). I'll update my answer with an description that vaguely fits your image, but be aware that it's best not to over-interpret PCA plots.

Thanks heaps Dr. Ryan, This helps. I appreciate your kind help. No above PCA is from targeted proteomics data of the sample.


text:
Which is the amino acids classification when analyzing the solvent accessibility of a protein?
I came across an article which has defined 3 categories: Buried [A, L, F, C, G, I, V, W], Exposed [ R, K, Q, E, N, D] and Intermediate [M, S, P, T, H, Y] when analyzing the solvent accessibility of a protein. But I could not find any reference validating the claim.
Another article has come up with 4 divisions:
B=completely buried (0-4% exposed)
b=partly buried (4-25% exposed)
e=partly exposed (25-50% exposed)
E=completely exposed (50+% exposed)
Which is authentic when dealing with solvent accessibility, and which should be used? Also for the second article [B,b,e,E], which amino acid corresponds to which class?

"All models are wrong but some are useful" is usually used in statistics, in this case, what do you want to solve? Can't you use both of them? If it is not explained in the 2nd article or in the supplementary information it will be hard to know by us to which amino acids are in each 4 categories. Did you looked there?

@Llopis yes I looked, no mention of anything. Since it wasn't mentioned. I thought it was probably a well known fact which is why they did not cite the source, hence I asked

There is some confusion in the question.
Solvent accessibility of a residue, i.e. the area accessible to the solvent (water), depends on the 3D structure of the protein. It's not a property of  an amino-acid (chemical compound such as alanine, leucine, ...).
There may be a correlation -- amino-acids with hydrophobic side chains are more likely to be buried, but in general amino-acids (chemical compounds) cannot be classified as buried or exposed.
I don't have access to the first article, but I tried running the prediction server described in the second linked article, and I got a result
GSGMQIFVKTFRDRLRNLLPQTITLEVEPSDTIENVKAKIQDKEGIPPDQQVLIFSRKRL
CCCCEEEEEECEEEEEECCCCEEEEEECCCCCHHHHHHHHHHHHCCCHHHEEEEECCEEC
EEEBeBbBebbEbEbEebEEeEbeBEBeeEeeBeeBbEeBeEeEEbEeEeBeBbbEEeeb

...

Porter (Secondary Structure):
H = Helix   (DSSP classes H, G and I)
E = Strand  (DSSP classes E and B)
C = Coil    (DSSP classes S, T and .)

PaleAle (Relative Solvent Accessibility):
B = very buried      (<=4% accessible)
b = somewhat buried  (>4% and <=25% accessible)
e = somewhat exposed (>25% and <=50% accessible)
E = very exposed     (>50% accessible)

In one place in the sequence you can see QQ, predicted as eB.
So clearly glutamine is NOT assigned to a particular RSA class.
Dividing RSA (relative solvent accessibility, 0-100%) into categories is arbitrary, I don't think there is any widely used convention. The easiest way would be to just use numeric values.


text:
estimate genome size: kmer-based approach from PacBio reads
Can anyone suggest a software/method for kmer analysis using PacBio reads (RSII)?
Something similar to Jellyfish, that I saw in a nice tutorial - but must be suitable for long, noisy reads. kmercounterexact from BBMapTools might also be an option, but again tricky with long, noisy reads.
The closest solution I could find was to estimate the genome size from overlap statistics (for the raw reads, the error corrected preassembled reads or by mapping the raw reads to the assembled contigs) by:

running an initial assembly using a small seed read length
then plot the pre-assembled read overlap histogram.

What do you think of this approach?
I would like to use this kmer analysis to estimate genome size, but for now I could not find much. I do not have companion Illumina reads - if so, I could try the approach proposed by Carvalho et al. 2016
I have reads & in-house genome assemblies (Canu) from different isolates of the same species of fungus to compare.

I don't think there is a method that would estimate a genome size using raw long reads.
The genome size estimates based on raw reads are done by fitting a model to kmer spectra (for instance Genomescope). The kmer spectra built from long reads are really messy due to the high error rate of long reads. That makes fitting of a model quite difficult. These methods assume small error rates. If the error rate is as high as 0.1 (which is normal for long reads), then the probability of having one kmer of 23 bases correct is ~8%. Meaning that if the coverage is 100x, just ~8x of homozygous kmers would be the real genomic kmers and it's very likely that kmer peaks of the sequencing errors, heterozygous kmers and homozygous kmers will just blend into one.
I can imagine, that if your coverage would crazy high (like 600x or more) the kmer spectra could separate errors and genomic kmers. You could also try to decrease kmer size to increase the proportion of correct kmers.
Another option would be to use a trick that they used to generate solid kmers for genome assembly using generalized de brujin graph using noisy reads. Maybe it would be better to take a look on the newer version of the assembler, that is called Flye.
But to be honest, I think for a long reads, a genome assembly is much more reliable way how to predict a genome size than any of the kmer tricks. I would simply assemble the genome, remap back reads, check the uniformity of coverage and conclude the genome size out it.
-- edit --
I tried to find some comprehensive evaluation of the assembly sizes compared to genome sizes, but I could not, so I just checked two examples, one big, one small and they both match quite well. Ammopiptanthus nanus was assembled to 823.74 Mb while kmer estimate was 889 Mb (ref). Plasmodium falciparum was assembled to 23.5 Mb while the original genome size estimate is 22.8 Mb. I know it's just two, but there is also a good reason to think that the genome size will be reasonable, feel free to add here more examples.
If you would be worried about contamination in the assembly, you can first run it though one of the methods to detect a contamination (like Blobology).

ok, I can get genome size from assembly, as discussed. thanks for the input

Using a genome assembly as a way to estimate the genome size is not reliable due to potential contamination causing assembly errors, possible failure to collapse some haplotypes, or assembly errors from inconsistent coverage.

@conchoecia I do agree, but the contamination problem is in kmer spectra analysis as well. But in my experience, long read assemblies I have seen so far were quite close to flow cytometry estimates and consistent with kmer spectra analysis. I will do some more research. Do you have any example where the assembly size was very biased?

@conchoecia better?

Hi @KamilSJaron, the only caveat that I had in mind is that both canu and FALCON do not handle haploid assemblies well if the animal is diploid and highly heteryzygous, like over 5% in the case of FALCON. For example, I am assembling ctenophore genomes which regularly have uncollapsed haplotypes due to their high heterozygosity. As a result, the genome assembly size is approximately twice what is expected based on a kmer spectrum histogram from Illumina reads. Most people probably won't be working with high-heterozygosity animals, but it is important to point out. That's all!

Ah, got you, I ve been there. I was assembling amphioxus genome that IS extremely heterozygous. I ended up setting the parameters so both haplotypes got assembled separately and it worked very well (assembled precisely the diploid genome size)

I have recently come across this submission from Wang et al. (2020), which is a method using the tool kmerfreq on error-corrected long reads (e.g. generated using Canu, FLAS, or LoRMA). It looks to be a promising approach.


text:
How to calculate Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation from protein (amino acid) sequences?
I want to code Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation in Python 3.6 from the scratch for amino acid sequences without using any package. Unfortunately, I have not been able to come across a proper literature which shows how these autocorrelation indices have to be calculated. Can anyone give a clear explanation as to how to calculate it or at least give a link which has got a clear description of how to calculate these?

Rajarshi Guha describes the formulas to calculate Moreau Broto autocorrelation, Moran coefficient and Geary coefficient and has a bibliography.
But why reinvent the wheel? It seems that the Autocorrelation module of the PyBioMed package does already what you want: it is a pure Python implementation without any external dependencies.

its for python 2.7, not 3.6


text:
Change sample names while using cummeRbund
I want to analyze RNA seq data using R package cummeRbund. Using this package I have a problem changing sample names. When I run this command in cummeRbund
samples(cuff)
sample_index sample_name sample_name parameter value
1            1          q1        <NA>      <NA>  <NA>
2            2          q2        <NA>      <NA>  <NA>
3            3          q3        <NA>      <NA>  <NA>
4            4          q4        <NA>      <NA>  <NA>
5            5          q5        <NA>      <NA>  <NA>

It shows q1,q2,q3,q4 and q5 as sample names for the five samples but I want to change the sample names and put following names for the samples Batis, Blue_silver, Chakwal_50, Local_white and Syn_22 for q1,q2,q3,q4 ad q5 respectively. To accomplish this I ran following command:
> samples(cuff) <- c("Batis", "Blue_silver","Chakwal_50","Local_white","Syn_22")
Error in samples(cuff) <- c("Batis", "Blue_silver", "Chakwal_50", "Local_white",  :
could not find function "samples<-"

but this command gave error as shown above.
Also the BAM files have same names for all the samples as mentioned above.
How could I change the sample names to desired names in R using cummeRbund?
Or I need to re-run cuffdiff and change the names somewhere else?

Don't know which file was used by cummeRbund for sample names but replacing q1,q2,q3 and q4  with Batis, Blue_silver, Chakwal_50 and Syn_22 for all the files in the directory containing cuffdiff output worked for me.
I used sed language for this as shown below:
sed -i 's/q4/Syn_22/g' *

Or more generally you can use
sed -i 's/sample_name_shown/sample_name_to_change/g' *

After that I deleted the sqlite object created by cummeRbund and again created this object using readCufflinks() in R.
This time sample names are changed as desired:
samples(cuff)
sample_index sample_name sample_name parameter value
1            1       Batis        <NA>      <NA>  <NA>
2            2 Blue_silver        <NA>      <NA>  <NA>
3            3  Chakwal_50        <NA>      <NA>  <NA>
4            4      Syn_22        <NA>      <NA>  <NA>


text:
To find p-value of compounds (pairwise)
I have next SNPs data.frame structure
> str(SNPs)
'data.frame':   1703 obs. of  4 variables:
 $ group: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ rs1  : Factor w/ 3 levels "D/D","I/D","I/I": 1 1 2 3 3 2 1 1 1 1 ...
 $ rs2  : Factor w/ 3 levels "a/a","a/b","b/b": 3 3 2 3 3 2 2 3 3 2 ...
 $ rs3  : Factor w/ 3 levels "G/G","G/T","T/T": 2 1 2 1 1 3 1 1 2 1 ...
 ...other rs

> head(SNPs)
  group rs1 rs2 rs3 ...other rs
1     A D/D b/b G/T
2     A D/D b/b G/G
3     A I/D a/b G/T
4     A I/I b/b G/G
5     A I/I b/b G/G
6     A I/D a/b T/T

For example, I noticed that rs5 4b/4a and rs6 G/G in group A very often occur together (see below). In group B they very often occur together too. So I want to know - is it statistical regularity or not.
I can create table with all pairs in both groups
> SNPs$rs5_rs6 <- paste(SNPs$rs5, SNPs$rs6)
> tmp <- table(SNPs$rs5_rs6, SNPs$group)
> tmp

              A   B
  4a/4a G/G   1  20
  4b/4a G/G  31  83
  4b/4a G/T  14  51
  4b/4a T/T   1   0
  4b/4b G/G  37 106
  4b/4b G/T  35 119
  4b/4b T/T  11  31

So, now I need compare (find p-value) group A and group B: 4a/4a G/G in group A vs 4a/4a G/G in group B, 4b/4a G/G in group A vs 4b/4a G/G in group B, 4b/4a G/T in group A vs 4b/4a G/T in group B and etc. Only this pairs/
How can I do that?
I need p-values for all
> tmp

              A   B
  4a/4a G/G   1  20 - p-value?
  4b/4a G/G  31  83 - p-value?
  4b/4a G/T  14  51 - p-value?
  4b/4a T/T   1   0 - p-value?
  4b/4b G/G  37 106 - p-value?
  4b/4b G/T  35 119 - p-value?
  4b/4b T/T  11  31 - p-value?

If you can restate your question in terms of explicit null and alternative hypotheses it will be much easier for people to help you determine an appropriate test to find p-values

@heathobrien chisq - it is my choice. Problem is not to choice test. Problem is to calculate it for pairs in tmp

If your question is how to apply a chisq test to the above table, the answer is chisq.test(tmp$A, tmp$B)

heathobrien, my question is how to apply chisq test to tmp and get 7 p-values:
compare4a/4a G/G in group A vs 4a/4a G/G in group B and get first p-value, compare 4b/4a G/G in group A vs 4b/4a G/G in group B and get the second p-value and etc.

OK. Now we're getting somewhere. See my answer

You can get standardised residuals from chisq.test and convert these to p-values, but it's not recommended. The frequencies of each genotype depend on the counts of all other genotypes, so the p-values aren't independent and they aren't corrected for multiple testing, but if you must:
pnorm(-abs(chisq.test(tmp[,c('A','B')])$stdres)))

hm... What other options are possible? whether the differences between the two groups were statistically significant for all pairs (rs5+rs6 and all pairs - 4a/4a G/G, 4b/4a G/G and etc)?

If I were you, I would just focus on the overall chisq p-value (corrected for the number of SNP pairs tested). The frequency of allele 2 depends on the frequency of allele 1 and the frequency of the genotypes depends on the allele frequencies, so I don't think it makes much sense to test each genotype independently

Note that all you're doing is looking at differences in linkage disequilibrium between groups. There are a lot of tools out there to that efficiently

chisq.test - not very good for values less then 5. I'v try to use fisher.test pnorm(-abs(fisher.test(tmp[,c('A','B')], simulate.p.value=TRUE)$stdres)), but I'v got the error - Error in abs(fisher.test(tmp[, c("A", "B")], simulate.p.value = TRUE)$stdres) :   non-numeric argument to mathematical function

o the p-values aren't independent and they aren't corrected for multiple testing, but if you must - I'm not must. I just think it is good solution to use chisq.test. I just need to know - is there any differencies in two grops

Fishers exact test doesn't give you residuals. Having values less than 5 is probably the least of your worries if you're trying to obtain p-values for each combination of genotypes.

In that case, I'm completely confused and do not understand what I can do with my data and how...

Let us continue this discussion in chat.


text:
How to get all chemical compounds that interact with a given protein target?
I want to get all Pubchem Chemical compounds (CID) that interact with a given UniProt target. I have a list of UniProt targets and I need to study the relationship between these targets and chemical compounds.
Thanks

The easiest solution to the problem of mapping chemical compounds to protein IDs is to use the ChEMBL compound ID to UniProt mapping.
Doing this with PubChem is not so straight forward in part because its downloadable data format is a hideous obfuscated and badly documented mess split over multiple compressed files which in turn contain multiple compressed files...... For a publicly funded project they should really do better.


text:
How is PCR duplication rate computed in scATAC-seq?
Reading Cusanovich et al. (2015) I encountered the sentence:

We mixed pairs of cell lines (HEK293T or HL-60 versus GM12878), performed combinatorial cellular indexing, and sequenced the resulting libraries to saturation (65% duplicate rate).

How is PCR duplication rate computed?

Cusanovich, D. A. et al. Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing. Science 348, 910–914 (2015).

The answer is in the supplementary methods page 3.

Finally, all fragments in the same library with duplicate start and
  end coordinates were removed using Picard


text:
Deciding which samples go in which batch
I have 370 samples to sequence, we probably will end up using only 96 samples per run (due to the barcode with the primers we'll use). This means running 4 batches. To minimize the batch effect I need to decide which samples are sequenced in which run.
I have several variables from these samples, some of them are categoric (sex, disease status, treatment stage, ...) and some of them are numeric (age, clinical score...).
To minimize the batch effect I wouldn't group all samples from the same category in one run, and I would distribute as uniformly as possible the samples according to the numeric variables. However, I could find many papers about adjusting, taking into account, but I couldn't find a paper about software to help designing the runs.
Is there any software to help designing each batch?
If there are several I would like to point out that I would define best as the one with easier user interface

I found the OSAT package which does almost what I want.

However it is specific for Illumina and for a limited number of plates: "plates that hold 2, 4, 8 IlluminaBeadChip chips and have 24, 48, 96 wells, respectively."
It doesn't take in consideration the numeric values:
library(OSAT)
inPath <- system.file('extdata', package='OSAT')
pheno <- read.table(file.path(inPath, 'samples.txt'),
                    header=TRUE, sep="\t", colClasses="factor")
pheno$AgeGrp <- runif(576, 0, 100)
gs <- setup.sample(pheno, optimal=c("SampleType", "Race", "AgeGrp"))
gc <- setup.container(IlluminaBeadChip96Plate, 6, batch='plates')
gSetup <- create.optimized.setup(sample=gs, container=gc, nSim=1000)

The minDiff package reported in Stats.SE.

It handles both numeric and categorical data, but it can only optimize for two nominal criteria.

I wrote experDesign which does everything I wanted

It can handle both numeric and categorical data for an arbitrary number of variables. It is not restricted to an even number of samples per batch.

I've had to go through this pain, and found the best way to do it was to simulate the model fit to start assessing the balance of covariates. I've got two notebooks available on Github, which should help you deal with this issue, although I've not thought about how to accurately assess continuous variables.
While it's not wrapped in a package, you can use these notebooks as a template to fit to your problem. I may wrap this in a package at some point, but I haven't had the chance to do so yet.
This Notebook will allow you to set up all your variables, max size of a batch (may vary due to multiplexing, library set up, plate size, etc). Bin them into appropriate batches, then create a design matrix, model fit, and test.
Apologies in advance that it's a bit of a mess, but if anything isn't clear, open an issue on the repo.

Interesting, thanks for sharing! However I don't understand why do you perform a test? If I understood the script it is for a mock data, so what is the purpose of this? To estimate what would be a completely random output?

That's a fair point, and you're correct, there's no reason for it. Generally, it was just for completeness so that nothing was widely unexpected when performing the test.

After trying it. I realized that this assumes even number of sample per each category (except for batch), which is not my case

Yeah, unfortunately you'll need to modify it to deal with uneven cases.


text:
Alignment using secondary and tertiary features of DNA
I m trying to align different sequences of length not greater than 50 bp. Could I incorporate additional information such as stacking energy, entropy, bonds etc as additional criteria in the alignment?
The sequences are very similar.

Welcome to the site. Did you some literature review? Which bonds are you referring to? Bonds between each pair of sequences? Or bonds within a sequence? I think it would be clearer if you explain what is your goal (and maybe it provides some key words to help find more information).


text:
STAR Indexing Diference for Small and Large Genome File Output
I have a quick doubt on the output of the Genome Indexing, I have used the STAR program along with genome .fasta file and GFF file.
Genome size is 3GB, here is the file output
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
genomeParameters.txt

I have another small Genome 60MB in size, I did the genome indexing, here is the file output
chrLength.txt
chrNameLength.txt
chrName.txt
chrStart.txt
exonGeTrInfo.tab
exonInfo.tab
geneInfo.tab
Genome
genomeParameters.txt
SA
SAindex
sjdbInfo.txt
sjdbList.fromGTF.out.tab
sjdbList.out.tab
transcriptInfo.tab

My point here is that, why I got the extra information for my small genome size, but I didn't get the same for the big size genome. I do apply the same procedure for the both.
here is the below information. Only difference I made for the large Genome size is (--sjdbOverhang 99 --genomeChrBinNbits 15) to reduce the memory, but the rest of things are same for small genome.
#!/bin/bash
NUM_THREADS=12
mkdir DB
STAR --runMode genomeGenerate --genomeDir DB \
    --runThreadN $NUM_THREADS \
    --genomeFastaFiles XL9_2.fa \
    --sjdbGTFfile XENLA_Frog.gtf+gff3 \
    --sjdbOverhang 99 \
    --genomeChrBinNbits 15

Could anyone give an idea, why there is different, I am new to this field, so I am wondering about the difference in this.
Thanks in advance.
Cheer
San

The files missing in the case of the larger genome are the actual index itself - it doesn't look like STAR completed. Do you have STDOUT/STDERR from the runs?

Here is my below syntax, I use it for It                                                     #!/bin/bash
module load star/2.5.0c
NUM_THREADS=22
mkdir Sangv92db
STAR --runMode genomeGenerate --genomeDir Sangv92db
--runThreadN $NUM_THREADS
--genomeFastaFiles /ifs/home/skm9/san/xlaevis/XL9_2.fa
--sjdbGTFfile /ifs/home/skm9/san/xlaevis/XENLA_9.2_Xenbase.gff3
--sjdbOverhang 99
--genomeChrBinNbits 15,                                                                                It run and gene the files out same like the above I mentioned.

here is the memory allocation I used for it, still not able to complete the indexing, ... qsub -pe threaded 4-20 -hard -l mem_free=400G -l h_vmem=560G -l mem_token=80G RunSTARIndex2.sh...

You didn't get extra information with the small genome, the index for the larger one just wasn't made. My guess is that you ran out of RAM. For mammalian sized genomes ensure you have at least 30GB of RAM, likely more like 60GB. Also, make sure to look for error messages.

qsub -pe threaded 4 -hard -l mem_free=60G -l h_vmem=80G -l mem_token=60G RunSTARIndex.sh...do you think this much is enough, although one node has 128GB....hope this time it will....me running now...Thanks a million in advance...

Thanks for the information, I do run the same but it seems, its terminating and giving the same output files. ..this is what I tried in the last ...qsub -pe threaded 4-20 -hard -l mem_free=100G -l h_vmem=120G -l mem_token=60G RunSTARIndex2.sh....Each Node -128 RAM, is their other possible way to Index it?

Some genomes require gobs of memory. If you point me to the files you're trying to index I can try and do that for you.

Thanks a million Devon Ryan, I tried all possible way to increase memory like ....qsub -pe threaded 4-20 -hard -l mem_free=400G -l h_vmem=560G -l mem_token=80G RunSTARIndex2.sh, but still fail, I don't see any error message, its seems like it terminate on the half way...Could you let me know your email, My be I will email the same to you.

Genome Indexing for Xenopus leaves, v9.2, here is the link for the database http://www.xenbase.org/other/static/ftpDatafiles.jsp, Please let me know, if you can  do the indexing for it, Let me know, how I can transfer the index file too, thanks a million in advance ...


text:
GDC data portal cancer subtype
I need to use Gene Expression quantification data (using HTSeq-FPKM-UQ) for breast cancer data on the gdc portal (https://portal.gdc.cancer.gov/repository) to classify breast cancer subtypes (luminal A, luminal B, HER2-like, basal-like).
However, since I don't know the portal, I'm having problems in associate the EXP data and the label.
I could get the label either by using the old TCGA-barcode and use a dictionary I already have from an article (but I don't know how to get the barcode from the portal / API) or by using clinical data (but I don't know how to find metadata that link the EXP and the clinical data).
Is there anyone more used to the portal who could help me?

Your question doesn't give enough details for even someone familiar with this portal to try and help, I think. Do you have an example of a dataset you want to interrogate, and how it is deficient in the details you require?

I personally find the GDC portal not very user-friendly. However, some of the data is available from other sources. For example, you can download TCGA-BRCA FPKMs and phenotypes from UCSC Xena here: https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Breast%20Cancer%20(BRCA)
All the sample info is already combined and the labels are the classic TCGA barcodes.

I would give a try with the R/Bioconductor package "TCGABiolinks" (https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html). The vignettes are quite helpful for various steps of data retrieval, including filtering based on user defined criteria.

http://gepia2.cancer-pku.cn/ try this

These data are in Experiment Hub in Bioconductor
library(curatedTCGAData)
brcaRNA <- curatedTCGAData("BRCA", "RNASeq2GeneNorm", dry.run = FALSE)

If you could expand your answer to explain how this helps to link data from the samples to each sample it would make a great answer


text:
Batch Effects in RNA Seq Sample
I am using the R (using EdgeR) for the RNA Seq analysis, I had few batch effect samples like Control vs treatment. Could anyone tell me the best way to remove the batch effects.
I have looked into the cluster dendrogram and MDS post (using EdgeR), the Control and Treatments are not grouping together, and their is differences in the different batch within Control and Treatment.
I tried this syntax to remove the batch effects for the rest of the analysis to get the DE Genelist. I feel it didn't able to remove the batch effects from the rest of the analysis.
# define the experiment
Batch <- factor(c(1,2,3,1,2,3))
Treatment <- factor(c(0,0,0,1,1,1))
design_noI <- model.matrix(~ 0 + Batch + Treatment)

# define our DGEList
San <- DGEList(counts=x,genes=gene_info)

A <- aveLogCPM(San)
y2 <- San[A>1,]
y2 <- calcNormFactors(y2)
logCPM <- cpm(y2, log=TRUE, prior.count=5)
plotMDS(logCPM)
logCPMR <- removeBatchEffect(y2, Clutch = Clutch)
plotMDS(logCPMR)

I can see the difference in MDS plot on removing the batch effects but I couldn't able to carry the same data flow into my rest of analysis.
# filter the lowly expresed genes and normalize the data
keep <- rowSums(cpm(San/logCPMR)>1) >= 3
keep2 <- removeBatchEffect(keep, Batch)
y_group <- San/logCPMR[keep2, , keep2.lib.size=FALSE]
y_group <- calcNormFactors(San/logCPMR)

I have mention the San/logCPMR, as San take as data without remove the batch effects, whereas I think logCPMR as batch effected removed data frame.
Any suggestion on this more useful to me.
Please let me know, any other information needed to deal with this batch effects.

I am not sure if I understand exactly what you're asking (and what you are doing in  your code). You seem to deviate from edgeR manual for paired or batch effect analysis. Your first 3 lines of code seem okay to me, but then you should fit your data to the model, including your design with batch effect. See 3.4 of edgeR manual. The removeBatchEffect results are only meant for clustering or visualization not the statistical analysis.

Thanks b.nota, My point is ...when their is batch effects, shown in MDS plot or  Cluster Dendrogram, re-visualise it by removing batch effects in MDS plot, But could we take the same batch effect removed data for further data analysis to get the DE genelist. As if we don't remove the batch effects from the data, that effects on the final DE genelist, and I get the less number of DE genelist.

No you cannot use removeBatchEffect results for DE analysis, follow the guidelines in manual and see answer of Devon.

The help page for removeBatchEffect() explicitly states that it should not be used in conjunction with differential expression testing, so don't do that. You already have the batch effect included in your model, so you're already correcting for it when testing for differential expression. If your MDS plots suggest that samples are variably affected by the batch effect then look into the combat() function from the sva package.


text:
How to efficiently extract coordinates from a `pybedtools.Interval`?
I'm using pybedtools to process data, and at some point, I loop over a BedTool object and use unpacking to extract coordinates from the Interval:
for interval in my_bedtool:
    (chrom, start, end, _, _, strand) = interval

Suppose I want to compute the interval length, as follows:
    length = end - start

This results in a TypeError: unsupported operand type(s) for -: 'str' and 'str'.
It happens that the unpacked end and start variables are strings.
Converting them to ints is probably suboptimal given that Intervals can be created by providing ints for the coordinates:
from pybedtool import Interval
interval = Interval("X", 0, 100, strand="+")

The start and end attributes of an Interval are ints.
I could do:
    (chrom, start, end, strand) = (
        interval.chrom, interval.start, interval.end, interval.strand)

But I think this looks a bit ugly.
Is there a better approach?

i do not know pybedtools but list(Interval("X", 0, 100, strand="+")) >>> ['X', '0', '100', '.', '.', '+']- shouldn't the conversion of integers to strings be considered a bug?

@Chris_Rands I was indeed surprised by this behaviour. As I say in my answer, I guess it might be related to how indexing works. It seems to rely on a string representation of the genomic interval: https://daler.github.io/pybedtools/intervals.html#indexing-into-interval-objects

I'm not sure, but the explanation of the behaviour could be found around these lines: https://github.com/daler/pybedtools/blob/master/pybedtools/cbedtools.pyx#L512

Or rather, this one: https://github.com/daler/pybedtools/blob/master/pybedtools/cbedtools.pyx#L256

Yes exactly, I was peaking at the same source code myself; so it's a design decision not a bug i guess, still seems odd to me (my guess is that keeping a homogenous list in terms of types might allow some performance benefits for cython)

The reason why unpacking an Interval results in strings instead of ints might be related to how indexing by position works.
Even without considering the string/int issue, it seems slightly more efficient to access the interval's attributes rather than unpacking them directly:
In [311]: %timeit [(chrom, start, end, strand) for (chrom, start, end, _, _, strand) in my_bedtool]
306 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [312]: %timeit [(interval.chrom, interval.start, interval.end, interval.strand) for interval in my_bedtool]
238 ms ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

To resolve the ugliness issue, attribute access can be made more compact through the use of operator.attrgetter:
In [313]: from operator import attrgetter

In [314]: get_interval_attrs = attrgetter("chrom", "start", "stop", "strand")

In [315]: %timeit [get_interval_attrs(interval) for interval in my_bedtool]
239 ms ± 26.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

FYI https://pastebin.com/3QBnVsKU


text:
A basic question about phylogenetic inference with Bayes-MCMC
I am new to the world world of Bayesian-MCMC inference and have been working my way through various texts. I am having some difficulty understanding how MCMC works. I understand that it is a way of sampling from the posterior distribution. What is puzzling me, however, is how the probabilities of the posterior distribution have been calculated. That is, if the Markov Chain is sampling from the posterior distribution, that posterior distribution must have probability values on some axis. How were those calculated?

In short, the posterior is a probability value between 0 and 1 which is calculated from Bayes rule, which In phylogenetic usage is given as the product of the prior probability and the likelihood divided by a normalising constant. The posterior probability allows you to state the probability of some hypothesis given the data at hand (as opposed to the likelihood parlance, which is probability of data given hypothesis).
In a phylogenetic setting, the sampler works (roughly) like this: propose some value of a parameter (topology, branch lengths, etc.) and use the specified prior and calculated likelihood for that parameter value given the data to calculate the unnormalised posterior. Next, propose a new value, and repeat the same step. Compare the two values, and if the posterior probability of the second is greater, accept this parameter value and repeat the cycle again. In this way the Markov Chain moves through posterior space until it reaches a reasonably stable optimal approximation of the posterior distribution. Here, we want to get as many samples as possible to maximise the accuracy of our posterior. The maximum a posteriori (MAP) value of the parameter and its credible intervals, etc, will come from the distribution of these samples.

Thank you so much! This is very clear and helpful. There is just one detail that I am hung up on, namely this step: "use the specified prior and calculated likelihood for that parameter value given the data to calculate the posterior." How can this be calculated without knowing the marginal probability/normalizing constant?

@Namenlos this is actually related to a property of the Metropolis-Hastings algorithm: because each step of the chain explores ratios of a value which is proportional to the posterior (that is, the likelihood multiplied by the prior; the numerator of Bayes' rule), the denominator of the posterior cancels out. Samples from the chain therefore give you an estimate of the posterior probability. I've edited my answer to reflect this.

Extra details can be found in this excellent introduction to Bayesian phylogenetics

Ah, I see. Thank you again. This is starting to make sense now.


text:
Omics data: How to interpret heatmap and dendrogram output?
How to interpret heat map and dendrogram output for biological data (omics) in words (when writing results and discussion)?
What should I consider (statistics behind?) and what is the best approach?
Here is one of my HM for proteomics data.
Script being used
col <- colorRampPalette(c("red","yellow","darkgreen"))(30)
png("HM1.png")
heatmap.2(as.matrix(df), Rowv = T, Colv = FALSE, dendrogram = "row", #scale = "row", col = col, density.info = "none", trace = "none", margins = c(7, 15) )
dev.off()

Heat Map output.
I would be greatful if you would give an example.

Why are you producing the heatmap in the first place? What question are you trying to answer?

Did you read the original paper on clustering by Eisen et al. in 1998? What have you tried already to understand more about why other people use heatmap clustering?

The dendrogram summarize the information of a group of values and sort them according to the similarity they have. It can be applied to both, samples and features.
The dendrogram allows to visualize features that are more similar together, usually revealing patterns that wouldn't have been seen otherwise.
In an article usually it is used something like: "we can see that X and Y are clustered together revealing that ... (they have something in common), while Z and A are clustered far"

You may be interested in reading up on heatmaps. For a history perspective (pre the biological introduction by Eisen et al. ) read The History of the Cluster Heat Map by Wilkinson and Friendly


text:
How to uninstall hmmpress from my conda environment?
I have installed hmmpress on my Mac OS X. It is here:
which hmmpress
/Users/user/anaconda/bin/hmmpress

It is available in my root environment:
hmmpress -h
# hmmpress :: prepare an HMM database for faster hmmscan searches
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: hmmpress [-options] <hmmfile>

Options:
  -h : show brief help on version and usage
  -f : force: overwrite any previous pressed files

I don't remember exactly how I installed it but I thought that I used conda because it is in the anaconda folder. However, conda search hmmpress returns no packages. I thought it could be part of the hmmer package. Indeed, the hmmer package exists in the bioconda channel. However, if I try conda uninstall hmmer, it says:
conda uninstall hmmer
Solving environment: failed

PackagesNotFoundError: The following packages are missing from the target environment:
  - hmmer

Does anyone know how I can remove hmmpress? I thought of rm -rf the /Users/user/anaconda/bin/hmmpress but it does not seem like a clean solution. There might be dependencies elsewhere that I could also get rid of.

That's normally part of the hmmer package (at least version 3 on OSX), so do a conda list | grep hmmer and see if maybe it's hmmer2 that you have installed.

This time conda uninstall hmmer in the root environment did the job. The sequence of commands that led me there is shown below. Thanks Devon Ryan.
conda list | grep hmmer
hmmer                     3.1b2                         3    bioconda

conda uninstall hmmer
Solving environment: done

## Package Plan ##

  environment location: /Users/user/anaconda

  removed specs:
    - hmmer

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    conda-4.5.9                |           py36_0         626 KB  conda-forge

The following packages will be REMOVED:

    barrnap: 0.9-0        bioconda
    busco:   3.0.2-py36_6 bioconda
    hmmer:   3.1b2-3      bioconda
    prokka:  1.12-3       bioconda

The following packages will be UPDATED:

    conda:   4.5.8-py36_1 conda-forge --> 4.5.9-py36_0      conda-forge
    openssl: 1.0.2o-0     conda-forge --> 1.0.2o-h470a237_1 conda-forge

Proceed ([y]/n)? y

Downloading and Extracting Packages
conda-4.5.9          |  626 KB | ################################################################################################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

Checking that hmmpress is not available anymore:
hmmpress
zsh: command not found: hmmpress


text:
Best practices for deciding if two structural variants are actually the same variant?
I know that I can use bcftools isec to compare two VCF files containing single nucleotide variants (SNVs); it will generate four possible output files: one with all of the unique called SNVs in file 1, one with all of the unique called SNVs in file 2, one with all of the SNVs from file 1 that were also found in file 2, and finally one with all of the SNVs from file 2 that were also found in file 1.
Is there a similar type of tool that will allow one to compare structural variants (SVs) between two VCF files?  I am about to tweak one of my Perl scripts (that currently uses the VCF.pm Perl module to parse VCF files) to perform this calculation (because thus far I have not found any tools designed for SVs) but I thought I would ask in this forum in case anyone else has already invented this wheel.

It looks like no one has answered yet.  Maybe others would find it easier to help you out if you provided a minimal example of the VCF files that you are trying to compare?

The tool that you want to use for this task is bedtools window.


text:
case/control phenotype prediction - population stratification with R
I'm struggling these last days on a case/control phenotype prediction problem based on genotype (individual x SNPs with more than 10k individuals and approx 18k SNPs encoded as 0,1,2).
Thanks to PCA analysis (with SNPRelate package), I observed that there are sub-populations in my training set. I suspect that there is a bias in my classifier due these sub-population.
My current methodology to improve my classifier accuracy is the following:

build one model for each population,
during prediction phase based on test dataset -> assign each individual to one of the population identified previously. Apply the
model built for this population.

My difficulties:

which packages and methods in R would you advise me to assign a
cluster to each of my individual ? Could I use PC eigenvectors to
achieve such clustering ?
which packages and methods in R could allow me to associate training genotype to cluster prior to apply the specific model ?

Thank in advance for your help.

Welcome to the site. You talk about your method, but we don't have any introduction about it. Are you testing with cross validation? Do you have some samples reserved for the final method? What is your goal? It seems to be to classify case and controls, but then you start adding sub-populations. If you want to classify regardless of the sub-population you will need to find shared features between populations, or first find the populations and then classify them between case and controls.

could you clarify what methods have you tried? At the moment I find a bit unclear what to answer

Hi, when the model is built (final step of my pipeline), yes, I'm using cross-validation. My goal is to predict a case/control phenotype. I'm thinking to build models for each of the sub-populations. Which methods do you recommend to find populations in R ?

What you consider a good method? Could you list the method you found/tried to avoid suggesting the same methods? How big is your dataset ? (some tools/methods perform better in bigger datasets than other methods)

Dataset size is 30k individuals by 20k SNPs. For the time being, I didn't yet find a method allowing me to identify sub-population, then to be able to predict sub-population on training sets. Any proposal will be welcome :)

@aure_bobo It's good to [edit] your questions with clarifications, so people don't have to dig in the comments to understand.


text:
The algorithm behind makeblastdb?

What is the algorithm/concept behind‚makeblastdb?

I know makeblastdb produces a .phr, .pin, and a .psq file, that represents an index or a data structure that is better optimised for searching a given query.

What is exactly stored in these files?

Seems like .phr and .pin are storing some metadata (based on the size) and .psq the actual index.
It can't be the table of all k-mers since the command doesn't take a value of k? (or maybe it is fixed to 3)
I could not find any direct discussion of the algorithm since google is flooded with people asking how to build such a database...
Maybe someone can give me a link or hint?

From looking at this publication "Database indexing for production MegaBLAST searches" which although deals with MegaBLAST and it's indexing method, it would appear that this was rolled into the NCBI C++ toolkit as the command makembindex so I'm inclined to think the current implementation must have been derived from this strategy. Looking at the relevant part of the current NCBI C++ toolkit I find no reference to makembindex but only makeblastdb which I assume has superseded the older command: blastdb Directory Reference
Ulimtately the true details of the index generation lie in the C++ source code for makeblastdb can be found here.  However following on from that I think the further details of the index generation may be in one of the following files which are includes for the main .cpp file:
seqdb.hpp
build_db.hpp
seqdbcommon.hpp
writedb.hpp
It would appear from seqdb.hpp that:
The index files have extension .pin or .nin and sequence files .psq or .nsq,
However it's not at all clear how the current implementation differs to the 2007 reference or where the .phr file is generated from!  Some degree of reverse engineering the C++ source code will be needed!

Thank you, i will read this! After a quick look, this datastructure is more or less a kmere list to position (hash/table), see 2.1.


text:
How to make year scale bar unit in figtree software?
I have generated phylogenetic trees for virus sequences by using beast software.
I observed in various research articles that the phylogentic trees also have a scale-bar (which mostly represents residue substitution per site). When I visualized the generated trees in FigTree v1.4.3, the scale bar unit (number) seems to be in thousands! While mostly the scale bar represents 0.1 or 0.5 such numbers.
I intend to add a scale bar to my year scale bar to phylogenetic tree, please advice me how to generate the optimal scale and value in accordance with the graphs.

Could you add a link to your current image to see what is some problem? But the problem is that you need to modify the scale of a bar or add a new one?

The scale bar is scaled to the range of the branch lengths in your tree file. If the branch lengths are in units of years but you want to show substitutions per site in the scale bar, you'll need to produce a tree file with branch lengths in units of substitutions per site.

In figtree time scale is showing in reverse order how to resolve this problem?https://ibb.co/j9WRke                          i attach link with this comments

In FigTree v1.4.3 there is some bug when i reverse time scale bar it didn't work....In  this version FigTree v.1.4.2 it working fine.


text:
Should PCA be standardized for gene expression?
This is a theory/good practice question more than a technical one. If samples are being plotted on a PCA projection of gene expression data, I'm wondering whether it is standard (and if so, why) to center and scale the PCs.
The reason I ask is that in this case the variables (genes) are being measured in the same scale, so I don't know if it would be necessary to center/scale data.

When the gene expression is scaled and centred you reduce the difference between genes.
Imagine you have one gene A that is highly expressed usually and has a standard deviation of 500 units compared to a gene B that is not much expressed and only have a standard deviation of 5.
In the scaled and centred genes both contribute the same because A usually is expressed 10000 and B is usually expressed 100 units. So, for both the standard deviation is 5% of their expression. Meaning that a variation in one is as important as in the other.
If not scaled (and centred) the first gene A will contribute more to the variation than gene B, because the expression variation in absolute numbers is bigger.
Both are used in publications, but I think the scaled and centred is more used, because the first dimension reflects "better" the differences between samples. Of course you can do both and select the one that is better to show what you want.

...on the other hand, if you scale and centre, a gene with expression 100 becomes as relevant as a gene with expression 10000 which may be undesirable. In the past, I followed the reasoning of the OP and opted against scaling and centring and I used log-transformed expression values (TPM or RPKM). Does this make sense?

@dariober The reasoning I use is, if a gene is usually expressed 1000, it needs this high expression to function, while the other doesn't need to be that high. So, why would I give more relevance to a gene because it needs more expression to keep the cell alive? What I want to see in a PCA is the similarities and differences, not if a gene is differentially expressed or not. About which units to use is (I think) another debate, as the normalization applied will change the apparent similarity between samples


text:
How to assess the quality of assembled .fasta genome files?
I have assembled 3 .fasta files from contigs infastq format of 3 different Homo sapiens.
I would like to see if the assembled genomes which are in .fasta format are in good quality or not. In addition, I would be happy to see heritages of genes from the maternal or the parental genomes.
More specifically, I have 3 .fasta files of a mother, father and daughter.

How can I see they are indeed the biological parents?
How can I study what genes were inherited from whom of the parents?
How can I assess the quality of the assembled genomes? Is there a standard way to do that?
Is there a simple frequency analysis of bases n-grams I can use?

Are these for humans? I imagine it'd be easier to work with VCF files than assembled fasta files.

@DevonRyan they are humans

What gene is this? Or it is the whole genome in a single fasta file? If the later case, if you explain the pipeline you use it might be easier to get an answer via VCF or tandem repeats.

@Llopis whole genome

To expand on the comment by Devon Ryan, what you really need are (A.) bam alignment files for each individual and (B.) VCF files for each individual showing all of their genomic variants (with respect to the reference human genome).  You mentioned that you had fastq files and that you converted them into fasta files.  Do you know if the bam files and VCF files are also available for download anywhere?
If the bam files are not available you will have generate them using sequence aligner software (for example: bwa-mem).  With the bam files you can calculate the average depth of coverage (and many other things).
If the VCF files are not available then you will have to generate them from the bam files using variant calling software (for example, the GATK package from the Broad Institute).
In terms of assessing the quality of the sequence runs, a commonly applied tool is name FastQC, but once again, it uses fastq formatted files as input, not fasta files.
The main use for fasta formatted files these days seems to be for compact storage of various reference genomes from different genome sequencing projects (using gzip or similar type of compression).
With the VCF files for each individual in hand you should be able to begin addressing your other questions regarding the frequency of shared variants, etc.

.bam and VCF aren't available. Is there any statistics that a valid genome should follow. Like frequencies of bases or ngram bases? Can you recommend of tools to generate VCF and .bam files from .fasta?


text:
Python API for working with ENSEMBL genomes
Dear bioinformaticians,
I could find two python APIs to mine ENSEMBL genome databases. There is a pyensembl which is more popular and then there is genomepy which is less popular. Following is my experience with both.
1. pyensembl:
pyensembl is a fast and powerful API to mine ENSEMBL genome databases.
I have been using it for a while and in my experience it works very well. The name of the package may make you believe that it is developed/maintained by ENSEMBL. That is not true. Still, being an open source package with a realistic potential to be the default API for mining ENSEMBL genome databases.
However, there are some issues related to it's development that users may find disconcerting.
For example, documentation. I could only find a comprehensive api documentation for release 0.8.10. Current release is 1.1. Not even release notes are available. Therefore users have to rely on few lines of documentation that is available in the README.md.
Also, as far as I know, it is not peer reviewed even though it is an excellent analysis resource.
2. genomepy
genomepy loads all of the genome files into the RAM, making it immediately unusable, especially if you are working with larger genome such as human genome.
So my question is: Is there any other better python API for working with ENSEMBL genomes?

"That makes you slightly worried about it's reliability."  You have the source code, you can see what it does and you can copy and modify as you wish. Why do you think it is unreliable? You can improve the documentation and suggest changes, but it looks like it is documented well with docstrings although the documentation has not been rebuild in a pdf format in a couple of years.

@Llopis I agree with both points. (1) I just wonder why it is not published. It is is a very useful resource. I will remove the reliabilitypart from my question. Thanks.
(2) Conventionally, for contributing (making pull requests), a potential contributor has to understand the code and to understand the code, the he/she has to go through the  documentation (preferably the latest one). So that's why I feel documentation is an important issue. May be I will build the documentation myself and contribute.

If we wrote papers on every tool we create we'd never have time to do anything productive. That and the limited benefit of actually submitting a paper on a tool/package is the reason that most never get published in that way.

The correct API for Ensembl is the Ensembl REST API which is updated and maintained by Ensembl, and language agnostic.

Hi @Emily_Ensembl, I wonder if there is a easy-to-learn kind of tutorial on Ensembl REST API. I want to use it to carry out general tasks such as (1) converting gene ids to symbol and other way around and (2) fetching a genome sequences of a particular species (assembly and ensembl release) within specific coordinates.

We've got an online course here.

Now there's another option, the ensembl_rest module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here.
To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.

This looks pretty neat and actually a lot more cleaner than the default python option in the EnsEMBL page. Great work!


text:
How to force pilon to correct a reference with low coverage mapped reads?
I have mapped some reads to a reference and want to change SNPs in the reference to SNPs found in the reads, even at low coverage. For example, here is a screenshot of one reference in IGV that I would like to correct.

When I use pilon to try to correct the reference, however, no bases are ever changed, even if using the option --mindepth 1.
Here is the call that I am using:
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --frags {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1

How can I use even low-coverage information to correct my reference?

Pilon will not consider fragments that are not proper pairs if the --frags option is used to pass the alignment .bam file. The reads mapped in the above IGV screenshot are unpaired.
Instead, if you want Pilon to consider all reads that map as legitimate and usable for correction, pass the .bam file with the --unpaired option instead.
To above corrected example would be:
java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --unpaired {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1


text:
Generating list of known fusions for Arriba from COSMIC
The RNA-Seq fusion detection program Arriba features the ability to white-list known gene fusions, specifically the authors recommend using COSMIC as a source of these here. Specifically they state:

A comprehensive list of known fusions can be obtained from CancerGeneCensus in the section titled "Complete Fusion Export". Depending on the gene annotation that is used to run arriba, some gene names need to be adjusted.

Emphasis is mine to highlight the issue here!
The state the format of the input file they require to be:

The file has two columns separated by a tab. Each line lists a pair of
  genes. The order of the genes is irrelevant. arriba searches for both
  genes as the 5' end and the 3' end of a fusion. Lines starting with a
  hash (#) are treated as comments. Optionally, the file can be
  gzip-compressed.

However the format of the download from COSMIC is somewhat different:
Sample ID       Sample name     Primary site    Site subtype 1  Site subtype 2  Site subtype 3  Primary histology       Histology subtype 1     Histology subtype 2     Histology subtype 3     Fusion ID       Translocation Name      Fusion Type     PUBMED_PMID     ID_STUDY
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1271    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Observed mRNA   16784981
WWWWWWW WWWWWWW thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1272    CCDC6{ENST00000263102}:r.1_535_RET{ENST00000355710}:r.2369_5659 Inferred Breakpoint     16784981
XXXXXXX XXXXXXX thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16784981
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS                              16728576
YYYYYYY YYYYYYY thyroid NS      NS      NS      carcinoma       papillary_carcinoma     follicular_variant      NS      1216    PAX8{NM_003466.2}:r.?_PPARG{NM_015869.2}:r.?    Inferred Breakpoint     16728576

Note: numeric samples IDs have been deliberately replaced with letters here.
Using this file, as is, will cause arriba to spurt out many warnings that the gene names in the first two columns are not recognised:
WARNING: unknown gene in known fusions list: XXXXXXX
From the looks of the rest of my arriba output the currency of gene IDs its using would appear to be gene symbols taken from the gencode.gft I'm supplying it.
I'm thus assuming the format arriba requires for a known gene list to be:
SYMB1   SYMB2
SYMB3   SYMB4

Where the record separator is a tab character \t
My question is how can I reprocess the COSMIC "Complete Fusion Export" .tsv file to the format required for arriba to act as a white-list?

I managed to solve this with a one-liner:
grep -v 'Sample ID' CosmicFusionExport.tsv | cut -f 12 | grep -v -e '^$' | perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' | sort -u > CosmicFusionList.tsv

The various parts of the piped command are as follows:

Remove header:grep -v 'Sample ID' CosmicFusionExport.tsv |
Extract 12th field (Translocation Name):cut -f 12 |
Remove blank entries:grep -v -e '^$' |
Extract gene symbols from 12th column:perl -lne '/^(\S+?)\{\S+\d+}\S+_(\S+?)\{/; print "$1\t$2"' |
Remove duplicate lines:sort -u > CosmicFusionList.tsv
Optionally the final two stages could be:sort -u | gzip -c > CosmicFusionList.tsv.gz to make use of file compression (for what it's worth on a 3.9 kB file)

Since the authors of arriba have neglected to document this necessary data munging step, I've posted it here in Q and A format in the hope that it will be useful to others, and other solutions can be proposed.


text:
Using pysam with cython: htslib/kstring.h not found
I'm trying to learn to use cython to speed up some code based on pysam. My issue is not strictly speaking about bioinformatics, but rather about building tools using a bioinformatics library. I hope this is still relevant in this site (I hesitated posting it as a github issue to pysam, alternatively, but I was not sure that this was a problem with pysam).
In my project's folder, the bam25prime/libbam25prime.pyx file has the following content:
# cython: language_level=3
from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

In my setup.py, the setup function is called with this argument:
ext_modules = cythonize("bam25prime/libbam25prime.pyx").
When I try to do the build (running python3.6 setup.py build_ext), cython seems to miss some include paths, because I get the following error:
[1/1] Cythonizing bam25prime/libbam25prime.pyx
running build_ext
building 'bam25prime.libbam25prime' extension
gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -O3 -march=native -fomit-frame-pointer -fPIC -I/home/bli/.local/lib/python3.6/site-packages/pysam -I${HOME}/include -I/home/bli/include/python3.6m -I/home/bli/.local/include -c bam25prime/libbam25prime.c -o build/temp.linux-x86_64-3.6/bam25prime/libbam25prime.o
bam25prime/libbam25prime.c:497:28: fatal error: htslib/kstring.h: No such file or directory

If I manually re-run the gcc command, adding -I/home/bli/.local/lib/python3.6/site-packages/pysam/include/htslib, the header file gets found.
How do I get cython to correctly determine the path to the necessary header files?
I guess another way to formulate the question would be the following: How is cython supposed to determine what include flags to use?
Or maybe there is a problem with my pysam installation?
If I remember well, I installed it using pip3.6 --user install pysam. I have a ~/.pydistutils.cfg file containing the following settings
[install]
optimize=1
[build_ext]
include_dirs=${HOME}/include
library_dirs=${HOME}/lib
rpath=${HOME}/lib
user=1

Any advice appreciated.

Have you tried setting C_INCLUDE_PATH? That's usually sufficient.

Make sure htslib is in your include path

Shouldn't cython infer this automatically when a python library is installed? Or is it because htslib is not in itself a python library, although distributed together with pysam?

Based on a suggestion by Devon Ryan, I searched how to set C_INCLUDE_PATH for cython, and found this issue. It refers to a (currently not working) feature present in cython's documentation, that should let libraries inform cython of the correct localisation for the header files in the module's setup:
ext_modules = cythonize(
    "bam25prime/libbam25prime.pyx",
    include_path=pysam_get_include())

As noted in the issue discussion, this seems to have no effect, but a workaround exists, consisting in using the Extension class to store the header localization information.
So the general answer is: Some python libraries (including pysam) provide a function that returns the paths to their include directories. This can be used to inform cython, but has to be done "manually".
My setup.py now looks as follows:
from setuptools import setup, find_packages
# If you have .pyx things to cythonize
from distutils.extension import Extension
from Cython.Build import cythonize
from pysam import get_include as pysam_get_include

name = "bam25prime"
__version__ = "0.1"

# https://github.com/cython/cython/blob/master/docs/src/reference/compilation.rst#configuring-the-c-build
extensions = [
    Extension(
        "bam25prime.libbam25prime", ["bam25prime/libbam25prime.pyx"],
        include_dirs=pysam_get_include()),
    ]

setup(
    name=name,
    version=__version__,
    # [description, author, email, licence, ...]
    packages=find_packages(),
    scripts=["bin/%s" % name],
    ext_modules = cythonize(extensions),
    # ext_modules = cythonize(
    #     "bam25prime/libbam25prime.pyx",
    #     include_path=pysam_get_include()),
    install_requires=["pysam", "pybedtools"])

Then, python3.6 setup.py build_ext succeeds.
And for whoever might be interested, here are the other files of the project:
bin/bam25prime is a link to bam25prime/bam25prime.py, which looks like this:
#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
"""
When executed, this module collapses reads in a bam file to their 5-prime end
(reducing them to length 1) and writes the result in bed format.
"""

import argparse
from pybedtools import BedTool, Interval
from pysam import AlignmentFile
from .libbam25prime import collapse_ali

# pybedtools.Interval | pysam.AlignedSegment
# chrom               | reference_name
# start               | reference_start
# end                 | reference_end
# bool(int(strand))   | is_reverse
#
# Interval(reference_name, reference_start, reference_end, strand=strand)

def collapse_and_sort(alis):
    """
    Collapse elements from *ali* to their 5 prime end.

    The resulting :class:`pybedtool.BedTool` is sorted.

    :param alis: An iterable of :class:`pysam.AlignedSegment`
    :rtype: :class:`pybedtool.BedTool`

    :class:`pybedtool.BedTool` is an iterable:
    It can be iterated over in a for loop, or produce an iterator
    using :func:`iter` on it
    (see https://stackoverflow.com/a/28353158/1878788).
    """
    return BedTool(
        "\n".join(map(collapse_ali, alis)),
        from_string=True).sort(stream=True)

def main():
    """Main function of the program."""
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        "-b", "--bamfile",
        required=True,
        help="Sorted and indexed bam file containing the aligned reads.")
    parser.add_argument(
        "-o", "--out_bed",
        required=True,
        help="Bed file in which the results should be written.")
    args = parser.parse_args()
    with AlignmentFile(args.bamfile, "rb") as bamfile:
        collapse_and_sort(bamfile.fetch()).saveas(args.out_bed)
    return 0

if __name__ == "__main__":
    exit(main())

bam25prime/libbam25prime.pyx:
#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
# cython: language_level=3
"""
This library contains a cythonized version of a function to collapse aligned
segments to their 5' end.
"""

from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

def main():
    """Example use case."""
    print("No examples to show here.")
    return 0

if __name__ == "__main__":
    import sys
    sys.exit(main())

(Any advice or suggestions welcome: I'd like to improve my python package-making skills.)

For those interested, to use cython with pybedtools I had to define the following Extension (where glob is glob.glob and /home/bli/src/bedtools2-2.27.1/ is the location of the source code of bedtools).
    Extension(
        "bam25prime.libcollapsebed", ["bam25prime/libcollapsebed.pyx"],
        include_dirs=glob("/home/bli/src/bedtools2-2.27.1/src/utils/*"),
        language="c++"),