Repeatmasker output table I'm working with an unassembled genome of cephalopod and I would to find transposable elements. My idea is to use repeatscout to generate a library and repeatmasker to find the elements. After creating a library with repeatscout,  I run Repeatmasker but it generates a table like this:

 file name: OB_100DEC.fa
sequences:           100
total length:  209466439 bp  (184235452 bp excl N/X-runs)
GC level:         35.24 %
bases masked:   73100024 bp ( 34.90 %)
==================================================
               number of      length   percentage
               elements*    occupied  of sequence
--------------------------------------------------
SINEs:                0            0 bp    0.00 %
      ALUs            0            0 bp    0.00 %
      MIRs            0            0 bp    0.00 %

LINEs:                0            0 bp    0.00 %
      LINE1           0            0 bp    0.00 %
      LINE2           0            0 bp    0.00 %
      L3/CR1          0            0 bp    0.00 %

LTR elements:         0            0 bp    0.00 %
      ERVL            0            0 bp    0.00 %
      ERVL-MaLRs      0            0 bp    0.00 %
      ERV_classI      0            0 bp    0.00 %
      ERV_classII     0            0 bp    0.00 %

DNA elements:         0            0 bp    0.00 %
     hAT-Charlie      0            0 bp    0.00 %
     TcMar-Tigger     0            0 bp    0.00 %

Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %

Small RNA:            0            0 bp    0.00 %

Satellites:           0            0 bp    0.00 %
Simple repeats:  290711     15358198 bp    7.33 %
Low complexity:       0            0 bp    0.00 %
==================================================

* most repeats fragmented by insertions or deletions
  have been counted as one element

The query species was assumed to be homo
RepeatMasker Combined Database: Dfam_Consensus-20181026

run with rmblastn version 2.6.0+
The query was compared to unclassified sequences in ".../OB_100DEC_repeats_filtered1.fa"

This was my command to run repeatmasker:

./RepeatMasker -s -lib /home/RepeatScout-1.0.5/OB_100DEC_repeats_filtered1.fa /home/Workdirectory/OB_100DEC.fa

The file OB_100DEC_repeats_filtered1.fa is the output of repeatscout and looks like this:

>R=3 (RR=4.  TRF=0.000 NSEG=0.000)
TAAGGCGGCGAGCTGGCAGAATCGTTAGCACGCCGGGCGAAATGCTTAGCGGTATTTCGTCTGTCTTTACGTTCTGAGTT
CAAATTCCGCCGAGGTCGACTTTGCCTTTCATCCTTTCGGGGTCGATAAAATAAGTACCAGTTGAGCACTGGGGTCGATG
TAATCGACTTACCCCCTCCCCCAAAATTTCTGGCCTTGTGCCTATATTAGAAACGATTATT
>R=4 (RR=5.  TRF=0.122 NSEG=0.226)
ACACACACACACACACACACACACACATATATATATATATACATATATACGACGGGCTTCTTTCAGTTTCCGTCTACCAA
ATCCACTCACAAGGCTTTGGTCGGCCCGAGGCTATAGTAGAAGACACTTGCCCAAGGTGCCACGCAGTGGGACTGAACCC
GGAACCATGTGGTTGGTAAGCAAGCTACTTACCACACAGCCACTCCTGCGCCTATATATAT
>R=6 (RR=7.  TRF=0.134 NSEG=0.247)
TTGTTTCAGTCATTTGACTGCGGCCATGCTGGAGCACCGCCTTTAGTCGAGCAAATCGACCCCAGGACTTATTCTTTGTA
AGCCTAGTACTTATTCTATCGGTCTCTTTTGCCGAACCGCTAAGTTACGGGGACGTAAACACACCAGCATCGGTTGTCAA
GCGATGTTGGGGGGACAAACACAGACACACAAACACACACACACACATACATATATATATATATATATATA
>R=8 (RR=9.  TRF=0.391 NSEG=0.000)`

and it is composed by all the TE (not classified) find in all the scaffolds of the genome.

Can anyone explain why there are no TE?
The fasta file that I used is only a part of the whole genome, it is about 200 mb and it is composed by the 100 greatest scaffolds of my genome.
Maybe I made a mistake with the use of libraries. Anyway I installed RepeatMasker with rmblastn version 2.6.0+ and RepeatMasker Combined Database: Dfam_Consensus-20181026.
thank you

As far as I can tell, you use the -lib option for species where there is no existing RepeatMaster database. Therefore, all it does is attempt to locate versions of the repetitive elements from the -lib file in the target genome file. However, since your lib file doesn't have any information on what type of element each of your hits is, it looks like RepeatMasker is treating all of them as "Unclassified".

Note the following from the repeatmasker.help file you should have downloaded along with RepeatMasker:

  The recommended format for IDs in a custom library is:

  >repeatname#class/subclass

  or simply

  >repeatname#class

  In this format, the data will be processed (overlapping repeats are
  merged etc), alternative output (.ace or .gff) can be created and an
  overview .tbl file will be created. Classes that will be displayed in
  the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything
  with 'RNA' in it, 'Simple_repeat', and 'Other' or 'Unknown' (the
  latter defaults when class is missing).

This makes sense. RepeatMasker cannot know if an arbitrary sequence you give it is a SINE or a LINE or anything else. It will therefore rely on the categories present in the fasta headers. Since your file has no categories, it identifies everything in your lib file as "Unknown", which is presumably what the "Unclassified" section in your output refers to:

Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %


Calculating celltype proportion changes between conditions in single-cell data I want to test for significance between differences in cell count between conditions.

I have 2 conditions (N, T), 3 donors (samples) each, and four cell types (A-D) as shown in the example below.

I believe this problem is analogous to the differential expression test performed with edgeR or DEseq.

Example data:

Count   Sample  Celltype    Condition
555 P1  D   T
3829    P1  C   T
945 P1  A   T
222 P1  B   T
795 P2  D   T
5464    P2  C   T
839 P2  A   T
438 P2  B   T
629 P3  D   T
6426    P3  C   T
647 P3  A   T
457 P3  B   T
8   P1  D   N
5122    P1  C   N
3204    P1  A   N
1   P1  B   N
20  P2  D   N
6111    P2  C   N
2524    P2  A   N
26  P2  B   N
66  P3  D   N
7490    P3  C   N
5272    P3  A   N
3   P3  B   N

This publication used a negative binomial regression model of the counts of cells at each time-point.

From the methods section: "For each cell type, we model the number of lineage-labelled cells detected in each analysed mouse as a random count variable using a negative binomial distribution. The frequency of detection is modelled by using the natural log of the total number of cells of that type profiled in a given mouse as an offset. The time point of each mouse (0, 30 or 60 days post tamoxifen) is provided as a covariate. The negative binomial model was fit using the R command ‘glm.nb’ from the ‘MASS’ package. The P value for the significance of the change in labelled fraction size between time-points was assessed using a likelihood-ratio test, computing using the R function ‘anova’."

The method is available here as an R script.

The input x is a dataframe with the following columns:

x$variable # celltype
x$value # cell count
x$batch # condition+donor
x$timepoint # condition (e.g. adult, fetal)

Montoro, D.T., Haber, A.L., Biton, M. et al. A revised airway epithelial hierarchy includes CFTR-expressing ionocytes. Nature 560, 319–324 (2018) doi:10.1038/s41586-018-0393-7


Highly variable elements? Which type of genomic elements in the UCSC browser can be used as an example of the highly labile = highly variable = low conservation regions?

Transposable elements (TEs) are the classic example of rapidly evolving regions in eukaryotic genomes. These are often used as a proxy for neutrally evolving sequences, which are by definition not conserved/constrained. The UCSC genome browser includes various repeat tracks (e.g. RepeatMasker) that includes TEs.

Within protein coding regions, 4-fold degenerate sites are often used as a proxy for neutrally evolving sequences. These sites are where any nucleotide at this position leads to the same amino acid and therefore protein (synonymous mutations).

Potentially even more variable sequences would be those under positive selection, but there is not a single class of element that falls into this category i think.


P-value correction when evaluating correlation between gene and miRNA expression First of all I apologize without the question is very basic, I am taking my first steps in bioinformatics.
Data information
We are evaluating the correlation (using the Pearson, Kendall or Spearman method) between gene expression and miRNA expression using the corAndPvalue function of WCGNA.
The resulting structure would be a DataFrame containing all combinations between each gene with each miRNAs, containing the following columns:
Gene     miRNA      Correlation P-value
Gen_1    miRNA_1    0,959       0.00311
Gen_1    miRNA_2    -0,039      0.1041
Gen_1    miRNA_3    -0,344      0.0021
Gen_2    miRNA_1    0,1333      0.00451
Gen_2    miRNA_2    0,877       0.07311
...

Question
Considering the huge number of correlation tests we are going to evaluate, we need to adjust the p-values to avoid correlations due to chance. Bonferroni does not seem to be the best solution, so we would use Benjamini-Hochberg method (BH). The question is:
The BH correction for the Gen_1 x miRNA_1 combination, should consider the p-values of all combinations that include Gen_1 (Option 1), or should consider all the p-values of all the genes x miRNA combinations (Option 2)?
For example, let's assume an expression dataset of 20,000 genes and another of 15,000 miRNAs
Option 1:
To adjust Gen_1 x miRNA_1 we would use 15,000 p-values (Gen_1 x miRNA_1, Gen_1 x miRNA_2, ..., Gen_1 x miRNA_15000).
Option 2:
To adjust Gen_1 x miRNA_1 we would use 300,000,000 p-values (Gen_1 x miRNA_1, Gen_1 x miRNA_2, ..., Gen_1 x miRNA_15000, Gen_2 x miRNA_1, Gen_2 x miRNA_2, ..., Gen_2 x miRNA_15000 and so on).
Suplementary question
Documentation of the method fdrcorrection from Python Statsmodels library suggests that for negative correlations (that could be frequent in a mRNA x miRNA correlation analysis) Benjamini-Yekutieli would work better; is that right? Or Benjamini-Hochberg method would be appropiated for this case?
Any kind of help would be much appreciated, thanks in advance!

I made the same question in CrossValidated forum and got an excellent answer!
The important part:

You need to correct for all of the comparisons you are doing. So if
that's 300,000,000 comparisons you need to correct for that many
multiple comparisons.

For more information check the answer in the link above


Documentation for Transfrag class codes (cuffcompare) Is there any extensive documentation or description for each class of Transfrag class codes as reported by the cuffcompare tool in the Cufflinks package?

Official doc might not be the best.

E.g. what does contained mean (class code: c)? Or, what is a generic exonic overlap (class code: o)?

You can get a bit more information by looking at the source code comments. Contained is for single exon transcripts that are primarily contained (>80% overlap) within an exon of an already annotated transcript. Generic exonic overlap is for cases where a transfrag exon overlaps that of an annotated transcript. The strand doesn't seem to matter in such cases.


can we show the orig.ident in separate windows by using Seurat I'm curious about can we show the orig.ident in sprat windows by using Seurat? Here are my code and the fig I got. Is it possible that I can separate each orig.ident into a single panel? DimPlot(whole,group.by = "orig.ident", reduction = "umap")

The DimPlot() function of the new version of Seurat, Seurat v3 has a split_by argument, which splits the plot based on the levels of the variable provided. You basically need to provide a column from the @metadata slot.

If you are using Seurat v2, you can have the DimPlot() function return the ggplot object created and then use facet_wratp() or facet_grid() to generate facets based on your variable of interest.


Find similar proteins with similar tertiary structure? I'm a mathematician studying knotted curves and am trying to find some applications of my work in proteins. Using https://www.rcsb.org/, how do I find proteins that have similar  (but perhaps not the same) tertiary structure? Preferably they would have similar primary structure and function...
Perhaps what I'm looking for is "similar proteins", if there is such a thing?

Also, is there a way in the database to select more "common proteins" (important or well studied, not just some rare exotic one (if there is such a thing)?

As a general question on the board, more detail is often good, because of the answers can be inaccurate due to misunderstanding of the biological context.

There are several ways to do this and at least 3 that I'm aware of.

Blast - sequence similarity searching the PDB database
Contact map identity, I use this approach, but its complicated - and its precision or utility is context specific, i.e. it could be loads of work and miss the vital detail you need.
Direct mapping one protein structure against another, this in particular is true for specific domains, or sites within the protein and the result is angstroms. Its not difficult and I assume this can be scaled to the whole protein.

Point 1 "Blackboxes"

By far the easiest is option one and I suspect I-Tasser (here) use this principle before instigating "threading". This package is considered "good" (others appear to consider it the best of its type), I believe it has done well in CASP tests (Critical Assessment of Techniques for Protein Structure Prediction). So the simple answer is dump your sequence in there an watch the output. Pyre is something similar. These are interesting tools but certainly not definitive, really there a starting point.

Point 1 "Formal Blasting PDB"

This is the easiest and most mature approach for your question, without knowing the details. I've run a test to ask the following question.

Which protein structure available on PBD is most similar to Japanese encephalitis virus (JEV) based on protein homology?

Step A
I've no specific interest in the virus BTW. The most robust Blast engine IMO is NCBI here . I dump a JEV sequence into the box, and click three important options (otherwise it doesn't work).

Database = protein databank (PDB), i.e. structural data
Organism "Japanese encephalitis virus (taxid:11072)": the taxid will come up automatically (if not just proceed as normal)
Next to the Organism box I click "Exclude", otherwise this whole approach doesn't work
Obviously I select "blastp", otherwise it definitely will not work

The output is/was here (but is now inactive). Summarising what it presented, E is so high its "0" (>2e-139) protein ID "77.59%" and just states "envelope protein" Chain A. I hit the link and the protein that comes up has a name which ends in "6A0P_A". This is the PDB reference and Percentage ID and PDB reference are the only two bits of info you need. The PDB reference needs the "_A" removing, thus its 6A0P. The _A just refers to chain A of a dimer.

Step B
You then go back to the Web site you mentioned and input the PDB ID you stated into the search bar. This results in this quite interesting page here. The structure is Usutu virus, which surprised me.

Point 2 . This is much harder and not only involves specific Python code, but a post unsupervised learning analysis and results in an average similarity of structure over the whole protein. There are a number of published algorithms for generating a contact map, or its easy to write your own. Its more advanced, its complicated and unless you've seen loads of them they not immediately understandable. Contact maps angle at a very specialist approach of predicting and understanding protein structure. It is not for the faint hearted and some biological details will be omitted, which may or may not be important for your purpose.

Point 3 . I don't have research grade expertise here, so I forget the exact button to click, but its quite easy via PyMol and other graphical viewers to perform spatial similarity for specific epitopes and this is something you need to learn within your sphere. Basically you will load two proteins into the image viewer and then they will be a "similarity aligning" function and you zoom in on your output to generate an answer in angstroms. You will need to learn this approach for the later parts of your project.

Summary, The orthodox way would be to perform point 1, select several targets and then examine them in detail via point 3. I definitely recommend direct Blast if you go down this route, because it is the most tractable and avoids potential unknown subsampling.  I personally would not recommend Point 2 in a general context because you have to consider an unsupervised learning approach towards processing a heap of contact maps. Ultimately, which ever approach you take you can never avoid Point 3 (which someone else does, in my case).

Note This may or may not be obvious to your but just because two sequence have the closest homology doesn't mean they will have the closest structural similarity.


Lifting from hg38 to hg37 I am having strange issues lifting this dataset to hg37 from hg38. I am using the HGDP WGS data from ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/. The data set split by chromosome, so I first converted them to BED then merged them to a single file using plink1.9. However, when trying to lift to hg37/19 using either liftOver  or CrossMap v0.5.1 I get strange errors. I tried both the chain from ensemble.org and from UCSC, I tried them both compressed and uncompressed. I am starting to think I need to try the conversion on the VCF files but these are so large I would prefer not to duplicate them unnecessarily.
Data
The data looks fine. If we look at the BIM:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail hgdp_hg38_chrpos.bim
Y   Y:56887228  0   56887228    C   T
Y   Y:56887285  0   56887285    C   A
Y   Y:56887316  0   56887316    A   G
Y   Y:56887463  0   56887463    A   C
Y   Y:56887491  0   56887491    A   C
Y   Y:56887503  0   56887503    T   G
Y   Y:56887583  0   56887583    G   A
Y   Y:56887631  0   56887631    T   C
Y   Y:56887837  0   56887837    A   G
Y   Y:56887844  0   56887844    C   T

The data are in chrpos variant names here but I had the same result using the original rsid.
Chains
Chains look fine, and I am able to read them into R without problems:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head GRCh38_to_GRCh37.chain
chain 1 1 248956422 + 10000 297968 1 249250621 + 10000 267719 2
167417 80249 50000
40302

chain 1 1 248956422 + 347968 501617 1 249250621 - 248779253 248932902 3
153649

chain 1 1 248956422 + 585988 12949384 1 249250621 + 521368 13009210 4
1044707 1 0
3709 3 0
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail GRCh38_to_GRCh37.chain

chain 1 chrY 57227415 + 56821509 57217415 Y 59373566 + 58967656 59363566 11183
85168 0 1
78580 0 1
63496 0 1
3158 1 0
30382 0 1
54147 0 1
80974

CrossMap
user@desktop:/media/luks8tb1/data/genomics/HGDP$ CrossMap.py bed GRCh38_to_GRCh37.chain hgdp_hg38_chrpos.bed hgdp_hg37_chrpos.bed
@ 2020-08-25 20:04:26: Read the chain file:  GRCh38_to_GRCh37.chain
Traceback (most recent call last):
  File "/home/user/.local/bin/CrossMap.py", line 166, in <module>
    crossmap_bed_file(mapTree, in_file, out_file)
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/mapbed.py", line 32, in crossmap_bed_file
    for line in ireader.reader(inbed):
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/ireader.py", line 27, in reader
    yield l.decode('utf8').strip().replace("\r", "")
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xdd in position 20: invalid continuation byte

liftOver
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.bed hg38ToHg19.over.chain hgdp_hg37.bed hgdp_hg37_errors.bed
Reading liftover chains
Mapping coordinates
Data format error: expecting at least 3 fields in BED file (hgdp_hg38.bed)

Trying the PED file with liftOver
I converted BED file to PED with plink1.9:
plink --bfile hgdp_hg38 --recode --out hgdp_hg38

Then:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.ped hg38ToHg19.over.chain hgdp_hg37.ped hgdp_hg37_errors.ped
Reading liftover chains
Mapping coordinates
invalid unsigned integer: "HGDP00001"

There's a few question on Google about this, but nothing that looks useful.
CrossMap does not support PED.

BED format can refer to the ".txt" file format that stores chromosomes and positions of genomic features as well as some other things like name. bed format can also refer to the plink format, which belongs to bed, bim, fam trio. The plink bed format stores all the genetic variants, but the program needs to refer to bim file for their positions.
Genome liftover obviously only needs the coordination of each SNPs, and it does not need individual genotypes at all. In fact, the bed format in PLINK does not provide any useful information. Knowing this means that providing the plink bed file is wrong. And this is pretty clear from the error message.
Moreover, since the three plink file always needs to be used together, most programs refer to them as plink file format. So whenever people say BED format, it is the txt file format.


WGCNA co-expression network analysis with less than 20 samples I have two cancer subtypes data. Subtype A is 14 samples and Subtype B is 23 samples. I'm interested in identifying the functions of some LncRNAs in the Subtype A group. For this I'm using all protein-coding genes and interested LncRNAs expression data (normalised counts).

From the WGCNA tutorial, I see that minimum number of samples needs to be 20. So, I guess using only Subtype A 14 samples for the analysis is not possible.

Can I use both Subtype A 14 samples and Subtype B 23 samples together for the analysis? In case if I use both, once I get the modules how do I know that the genes in modules are related to Subtype A or Subtype B? The pathways I get from the interested module is related to Subtype a or Subtype B? How this can be identified?

It's hard to say without knowing how different the subtypes are. If you do a common WGCNA, you may find modules related to the differences between A and B as well as to drivers of expression variation in A and B (common and separate). Since A has fewer samples, the variation will have to be stronger to be reflected by WGCNA modules.

You can run the WGCNA, find modules that vary the most in A samples (these may or may not also vary in B samples), then see whether LncRNAs belong to these modules and have a high kME (module membership).


How do I get the compartment Id and volume from SBML using libSBML Python? I searched the Python documentation for libSBML but couldn't find a staightforward answer to this question. For a given SBML model, how can I use Python libSBML to retrieve the values of the id and volume attributes for a given compartment?

There aren't any examples in the python docs on how to do this. The best bet is to search the C/C++ docs which are more extensive. There is a python example in the docs which shows how to retrieve the number of compartments in a model but it doesn't go further. The code below shows how you can retrieve the information requested, assuming you have a model called testmodel.xml:

import libsbml

document = libsbml.readSBML("testmodel.xml")

if document.getNumErrors() > 0:
   print("Encountered the following SBML errors:" + "\n")
   document.printErrors()
else:
   model = document.getModel()

   nCompartments = model.getNumCompartments()

   print ("Number of compartments: ", nCompartments)
   for i in range (nCompartments):
       comp = model.getCompartment(i)

       print ("Id = ", comp.getId())
       print ("Volume = ", comp.getVolume())


Pointing Cromwell to local EBS mounted storage on AWS Batch I'm trying to run the WDL tasks from the GTEx RNA-Seq pipeline with Cromwell using AWS Batch as a backend.
I store the STAR alignment index on an Elastic Block Store since my Cromwell server instance lacks sufficient disk space. The EBS Block is mounted to my Cromwell server at /data. The STAR index is located at /data/genomic-references/rnaseq/star_index_oh50.
Try as I might, I cannot get Cromwell to run the STAR task correctly.
My JSON inputs file:
 {
   "rnaseq_pipeline_fastq_workflow.prefix": "s3://mybucket/test/CL100148068",
   "rnaseq_pipeline_fastq_workflow.star.fastq1": "s3://anotherbucket/genotypedata/sample1.fq.gz",
   "rnaseq_pipeline_fastq_workflow.star.star_index": "/data/genomic-references/rnaseq/star_index_oh50/",
   "rnaseq_pipeline_fastq_workflow.star.memory": "32",
   "rnaseq_pipeline_fastq_workflow.star.num_threads": "4",
   "rnaseq_pipeline_fastq_workflow.star.disk_space": "200",
   "rnaseq_pipeline_fastq_workflow.star.num_preempt": "0"
 }

Error message:
tar: /data/genomic-references/rnaseq/star_index_oh50/: Cannot open: No such file or directory

Note that /data/genomic-references/rnaseq/star_index_oh50/ does exist and is readable.
The relevant command in the GTEx container is here:
tar -xvvf ${star_index} -C star_index --strip-components=1

If the files are stored locally, then why does tar fail to wrap the directory?

Somebody helped me with this one through other channels, but I wanted to make it publicly available:
Make a tarball of the STAR index. The option --strip-components=1 strips the first level of the path from the untarred items; thus, the expected tarball was made with something like
tar -cvf ${star_index} the_star_index_directory/

Then point Cromwell/WDL to the tarball, e.g. on S3. The GTEx container will unwrap the tarball on the fly during execution.


Duplicate genes with RSEM counts: Which one to choose? I have Ensembl ids in the first column and samples with RSEM counts data in other columns. I converted Ensembl ids to gene symbols. Now I see there are three genes repeated twice.

Gene   S1       S2       S3
COG8  804.07  1475.16  323.80
COG8  249.23   867.94  113.64
SDHD 1464.84  3508.02    75
SDHD  4.16     3.97       0
SOGA3 81.65   30.76    32.69
SOGA3 0.00    10.77    25.65

These are the Ensembl ids:

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

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

There is no one-to-one mapping of gene ids from one database to the other. Ensembl (who maintain Ensembl ENSG IDs), ncbi (who maintain EntrezGene IDs and RefSeq transcript ids) and HUGO (who maintain gene symbols/names) all have different ideas about what a gene is, which part of the genome belongs to which gene etc.

In fact, it's worse than that because HUGO gene IDs refer to the abstract idea of a gene, rather than a particular part of any chromosome.

My suggestion would be to leave the IDs as Ensembl IDs (as these actually refer to a physical part of a chromosome, which is what you are quantifying), right up until you have to report the results to someone who really NEEDS gene names. If you still have duplicates, you could sum, or pick the most interesting, but I would suggest just leaving the duplicates in there, and treat it as thought there are two genes with the same name.


Inferring a phylogenetic tree from BLASTn I am trying to infer a phylogenetic tree from a Blastn output and from what I have understood, what I should do is 1) extract the alignments and re-align them using Muscle, then 2) feed the .aln file into PhyML.
Is this correct?

Or should I be extracting the hsps from each alignment as suggested here https://biopython.org/wiki/Phylo? I don't quite understand why the hsps are extracted here, shouldn't the whole sequences be compared in the multiple alignment?

If you look into the get_seqrecs() function, you will see that this function is not extracting hsps to use in the following steps but the hsps are used in order to filter out low score BLAST alignments. The function accepts a blast alignment result and a threshold, and if the hsp score is above the supplied threshold the whole sequence (not just the hsp) would be the output. Your intuition is correct, the whole sequence should go in the MSA step.


Efficiently aligning a lot of reads on the same small reference sequence The context: I have a DNA-sequence coding for a protein, about 1500 bp in length. Using NGS, a lot of reads of (mutants of) this same sequence were acquired. All of these reads need to be aligned to the reference. We're talking about a lot of reads (100,000 - 1,000,000) of not too short length (120 - 300 bp). All of these reads can belong to different mutants of the template sequence, so the alignment is necessary to determine the exact sequence of every single mutant.

Currently, I'm just using Smith-Waterman-based local alignment to align every single read to the reference one by one. Yet I can't help but feel there might be a more computationally/time-efficient solution to this specific problem.

Maybe there exists an algorithm that's not very efficient for most alignment problems, but that becomes very worthwile if it has to map a ton of reads to the same place over and over again. For example, it might do some time-consuming operations on the short template that make it fast to align reads to it, but these operations aren't worth it if only aligning a couple of reads. That's just an idea, I don't know all the different techniques that are out there.

So, to recap: I want to align a lot (100,000 - 1,000,000) of 120-300 bp reads to the same short 1500 bp reference sequence. If anyone has any suggestions about an algorithm, or just a specific workflow, that is particularly suited to do this, it would be appreciated. I work in R so I can implement some things myself, it doesn't have to be a ready-to-use software package or anything like that. Thanks in advance!

I would suggest using a whole genome aligner like BWA or BOWTIE2 and swapping out the whole genome reference with your specific reference.

These alignment algorithms are designed to align large numbers of NGS reads against a single reference.  They can be used with any reference including short ones like your example or small viruses.


Correlate DEGs from DESeq2, EdgeR and Limma results I have a lists of DEGs identified by DESeq2, EdgeR and Limma.
I would like to correlate the the gene rankings in the lists to decide on a package to use in downstream analysis.
I am havig a few problems with this, would anyone be able to help?

Thanks

I would use the p-value/FDR which each method returns to rank the gene list for that method in order from 'most likely to be DE' to least likely.  You would then have three ranked lists of genes - you would put the actual rankings in for each gene - 1 for the top ranked, 2 for the next, etc.  You could then carry out, say Spearman's rank correlation on the lists of rankings to assess agreement between the methods.

You could also decide to keep all 3 methods, and accept genes as DE if, for example, at least 2 out of 3 methods give an acceptable p-value/FDR.

From my experience, correlation should be reasonably high at the top of the list, where the methods tend to agree more on the rankings, and when the correlation gets lower, further down the list, would be where I would start to question whether the genes are really differentially expressed.


tool for the prediction of impact of mutations on protein localization I need to know if the mutations present in the "position of interest" can have an impact on the protein localization ability. Is there any tool to make such prediction that can make my life easier?

I know the position of nuclear localization sequence in a protein and I also know it has a mutation, but I don't know how to predict the impact of this mutation on the localization of this protein ie if this mutation hindering protein's entry to the nucleus or not? Should an analysis of protein stability be useful it this case, as an alternative?

I am not interested in localization itself but the impact of mutations on the localization.

Do you have the protein sequence after the mutation available with you?
If so then maybe you can BLAST them and find out where the best hits are localized using GO annotation (cellular components). I am aware that these won't be your exact protein sequence but they can at least give you a rough idea.


Why does RATT create a directory with chromosome sequences and how can I stop it? I have run RATT to transfer the annotation from a reference strain to another strain. For this, I have run:

$RATT_HOME/start.ratt.sh embl ref.fa output/annotation/ Strain

This creates several things in the working directory that I do not wish:

a Sequences folder with the sequence in FASTA format of each chromosome
a copy of the reference sequence in FASTA format (followed by a mysterious number)
a copy of the query sequence in FASTA format followed by a mysterious number

I could not find any options in the RATT help to disable the creation of these files. Also, does anybody know why these files are created?

RATT seems to be a very simple wrapper script that unfortunately does not clean up its temporary files properly. Short of modifying the code yourself, your only option is to perform the cleanup manually after running RATT.

The temporary copies are made because (a) the query is not an exact copy: RATT performs some cleaning. And (b) it needs to extract a FASTA file from the EMBL input for use in downstream tools.

The “mysterious” number behind the file names is simply the process ID of the RATT process. This is appended to make these intermediate files more or less unique (presumably to prevent overwriting existing files).


Why don't all reads have adaptors? I got NGS reads back from sequencing platform. I check for adaptors and trimmed them. But I realized only a fraction (eg 30%) have adaptors... why not all of them?

Thanks

Things like this might depend on your specific library prep, but in general: sequencing starts at the end of the adapter, not before it. You will only see adapters if you sequence through the entire fragment into the adapters on the other side. If your fragment is long enough you won't see any adapter, and of course, that's the desirable outcome of your library prep.


How to subset genes and its nested features from a GFF file using a gene list I would like to subset a GFF file (gene and nested features) from a gene list.
The GFF file looks like this
##gff-version 3
Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903

And a have the target genes in a list such as
Genecv11033552m
Genecv11003131m
Genecv11036683m
Genecv11012576m
Genecv11003654m
Genecv11012587m

I know that is possible to subset the gff using grep -f gene_list.txt <gff_file>. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five_prime_UTR, CDS, three_prime_UTR).
This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.
Any ideas?

Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython.
Dealing with nested data in bash or with some fancy awk script is probably going to trickier to write than it's worth.
You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with mRNA in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.
If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.
Repeat until you reach the end of the GFF file or your list of target genes is empty.


Rule Extraction from nnet results I used a script in R language that uses nnet library to predict promoter bacteria and i would like to know how to extract rules from this neural network results.

As my input of the neural network i have n positive examples of promoter and n false examples. I use as input these examples and the FASTA file with the genome of the bacteria. As a result, i have a value of 0 to 1 for each network tested corresponding to its learning. I want to discover roles that can improve my network in future experiments in the same bacteria.

Which algorithms or softwares could i use?

R's nnet package only supports fully connected neural networks with one hidden layer. This is the most primitive type of network. I doubt it will work well for promotors finding. In addition, such network won't give you useful interpretations.

If you want to explore neural networks, you should use a one-dimension convolution layer as is described in this paper. This layer effectively represents a positional weight matrix (PWM). You can know which motifs are used. With backtracking, you can also identify the coordinates of motifs. To deploy such a model, you need to be fairly familiar with deep learning and to learn one deep learning framework. For framework, you may start with keras. Once you get used to keras, you can implement your full promotor finder in less than 100 lines of python code. Alternatively, you may try dragonn. It is supposed to simplify deploying models for DNA sequences. I have no experience with it, though.

For promotor finding, it is also worth trying traditional methods. It is interesting that few/no neuralNet-based publications have evaluated traditional methods, probably because many ML people know little about classical motif finding. These methods could work well, too.


How to count reads in bam per bed interval with bedtools I recently installed Ubuntu 16.04 (because I was still using 12.04). But it seems my bedtools scripts don't work properly anymore. I can't figure out how to use the new bedtools for my old ways. What I want to do is get the number of reads from a bam file, per interval from a bed file. It was very simple with my old version of Ubuntu and bedtools:

bedtools coverage -abam file.bam -b All_peaks.bed > file.cov.txt

But this doesn't seem to work anymore. I used to get small files as results (~3MB), but now huge files are created (~4GB) including read names??

I am using bedtools v2.25 now.

The order of -a and -b switched at some point. You want:

bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt

For reference, this is the end of the help output in version 2.25:

Default Output:
     After each entry in A, reports:
       1) The number of features in B that overlapped the A interval.
       2) The number of bases in A that had non-zero coverage.
       3) The length of the entry in A.
       4) The fraction of bases in A that had non-zero coverage.

And this is the equivalent output from version 2.19:

Default Output:
     After each entry in B, reports:
       1) The number of features in A that overlapped the B interval.
       2) The number of bases in B that had non-zero coverage.
       3) The length of the entry in B.
       4) The fraction of bases in B that had non-zero coverage.

Update: The change in behavior happened in version 2.24:

  We have changed the behavior of the coverage tool such that it is consistent with the other tools. Specifically, coverage is now
  computed for the intervals in the A file based on the overlaps with the B file, rather than vice versa.


Sort reads in BAM file based on presence of specific deletion? I have an indexed BAM file containing long-read sequencing data and I'd like to split the reads contained within into those with a known deletion and those without the deletion (I have the deletion coordinates available to me) when mapped against the hg38 genome. Ideally I'd like to end up with two BAM files (one with reads containing the deletion and one without) for downstream analysis. I've searched a bit online and haven't really found any tools that do exactly that, but this seems like a common enough task that there has to be something.
Thanks in advance for any help/guidance!
UPDATE
In addition to @gringer's great answer, I decided to write a small python script myself that solves my issue using pysam, the fruits of which can be found here. If anyone has any feedback or suggestions on it I'll happily take them, as I'm still very new to this. It's undergone very little testing, so user be warned.

The CIGAR string should tell you if a read has a deletion in it.  So filter the bam files based on that.


Retrieve ID ligand from PDB file I have thousands of PDB files and I need to extract the following information from each of them: an ID or the name of the small molecule (ligand) present in the complex.

Is there a way to do that with BioPython or similar tools? Taking this PDB entry as an example, the tools should return 7FM.

I found out that the package PyPDB has the function get_ligands that does exactly what I'm looking for.