text:
How to iterate protein sequences using amino acids?
I was working on certain program using python and I have been using a protein aligned sequence file in two formats, phylip (.phy) and clustal (.aln).
Example clustal file:
CLUSTAL 2.1 multiple sequence alignment

    Homo            -------------------------------------------MEMQDLTSPHSRLSGSS
    Macaca          -------------------------------------------MEMQDLTSPHSRLSGSS
    Canis           REWSSARPERSKGRRKPVDAAAVSAVQTSQTSSDVAVSSSCRSMEMQDLTSPHSRLSGSS
    Mus             -------------------------------------------MEMQDLTSPHSRLSGSS
    Rattus          -------------------------------------------MEMQDLTSPHSRLSGSS
                                                               *****************

Example phylip file:
     3   120
Apis       MPWIKCLGGG GGKARRGKPL SVSQPIHLLV KSEFEPQCHV FRTKQLRKPR
Branchiost ---------- ---------- ---------- ---------- ----------
Nematostel ---------- ---------- ---------- ---------- ----------

           PCHLCHQAVI KQASCCRVCK YICHKTCEDK GARFRDTSQQ IHKQWQADPA
           ---------- ---------- ---------- ---------- ----------
           ---------- ---------- ---------- ---------- ----------

I want to navigate through the file in a way that ;

If the columns of the file have same amino acid, it saves and prints the amino acid that exists in the column.
I need to know how can I search the columns by using the amino acids, if particular group of amino acids exist in that column or not. For example, if 'STA'or 'HY' or 'FVLIM' or 'NDEQHK' (with no '-') exists in any of the columns, then put a '@' in the file at the end of that column.
Following is the code I have been trying to manipulate the file with the following program:

Code below
alignment = open ("example.aln")
a = b = c = 0
for col in (alignment):
    num = len(set(col))
        if num == 1 and col[0] != '-':
            print (num)
        a += 1
        elif num > 1 and '-' not in col:
        print (num)
        b += 1 elif '-' in col:  # assumes 1 or more dashes
        c += 1
        print('Number of columns with the same amino acid: {}\n'
        'Number of columns with at least 2 amino acids (no gaps): {}\n'
        'Number of columns with one/more gaps: {}'
        .format(a, b, c))

These variables return the occurrences of amino acids in the columns of the file but I don't understand how to search using the amino acids.
What is the best way to iterate over the file and write into a new file with results using python? Should I use python switch statement?

Hi Sidra, as it was the block of code was unreadable since it seemed like many of the lines had been merged together. Could you edit your post to make sure that the code displayed above exactly matches your program? It would also help if you make your code a working minimal example where you provide a small testcase "example.aln". https://meta.stackexchange.com/questions/22186/how-do-i-format-my-code-blocks

Cross posted at Biostars: https://www.biostars.org/p/336237/#337089

In response to question 1, save the following as a python3 script (e.g., aln_parser.py):
#!/usr/bin/env python3

from Bio import AlignIO
import sys

aln = sys.argv[1]

algnmnt = AlignIO.read(open(aln, 'r'), 'clustal')

for col in range(0, algnmnt.get_alignment_length()):
    proteins = set(algnmnt[:,col])
    if len(proteins) == 1:
        print(f"Postion {col}: protein {''.join(proteins)}")

Run the script to query your alignment (clustal.aln).
python3 aln_parser.py clustal.aln
To save the output to file:
python3 aln_parser.py clustal.aln > results.txt
Requires python3 and Biopython (pip3 install biopython).

Looks good. You could also add a hash-bang #!/bin/env python3 to make it executable.

question 2 has still been the problem, if anyone can help.!

@SidraYounas I would try to help, but I don't fully understand the question.

@user3479780 You can see the code sample above, that I have used. In this, I am iterating over the columns of the files. All I was trying to use the way if I can search for the presence of some particular amino acids for each column.

@user3479780 I have been trying to delete these columns which have a single amino acid, but I am not being successful. Can you help?


text:
Genomic coordinates of predicted binding pocket
I would like to know the genomic coordinates of a proteins ligand-binding pocket.  For example, on this PDB page you can see that there is xray crystallography data for how CYP2C9 interacts with two ligands. I would like to know which amino acids form the space on the protein that binds the ligands, then translate (de-translate?) that back to the nucleotide sequence. Is anyone familiar with either a database that has this information or how to extract the ligand interaction data from the PDB files?
Thank you

It turns out there is a tool that does this beautifully.  I used the proteinToGenome function from the ensembldb package.
The code looks like this:
# Define the Ensemble database object for the
edbX <- filter(EnsDb.Hsapiens.v75, filter = ~ seq_name == "X")
# Set the first and last amino acid you want to convert and the ensembl protein id
GAGE10_prt <- IRanges(start = 5, end = 9, names = "ENSP00000385415")
# Perform the conversion.  GAGE10_prt is an object that has the coordinates
GAGE10_prt <- proteinToGenome(GAGE10_prt, edbX)

That PDB page includes a publication reference and abstract, containing the following text:

Mammalian CYP450s recognize and metabolize diverse xenobiotics such as
drug molecules, environmental compounds and pollutants.
...
Here we describe the crystal structure of a human CYP450, CYP2C9, both unliganded and in complex with the anti-coagulant drug warfarin.

Here's the genomic context for CYP2C9, with Ensembl page here.
The binding pocket images from the paper here and here indicate that the following amino acids are potentially involved in the interaction:

Arg97
Phe100
Ala103
Val113
Phe114
Asn217
Pro367
Phe476

Using the filtering from the UCSC table browser (genome: Human; group: Genes & Gene Predictions; Basic table; Region: genome; Identifiers: CYP2C9; Output format: CDS fasta alignment), I get the following amino acid sequences respectively for the exons in CYP2C9. I have added in the protein reference range to the end, with the help of the R cumsum function applied to the region lengths, cumsum(c(56,54,50,54,59,47,63,47,61)):
>ENST00000260682.7_hg38_1_9 56 0 0 chr10:94938683-94938850+ [1-56]
MDSLVVLVLCLSCLLLLSLWRQSSGRGKLPPGPTPLPVIGNILQIGIKDISKSLTN
>ENST00000260682.7_hg38_2_9 54 0 1 chr10:94941858-94942020+ [57-110]
LSKVYGPVFTLYFGLKPIVVLHGYEAVKEALIDLGEEFSGRGIFPLAERANRGF
>ENST00000260682.7_hg38_3_9 50 1 1 chr10:94942192-94942341+ [111-160]
GIVFSNGKKWKEIRRFSLMTLRNFGMGKRSIEDRVQEEARCLVEELRKTK
>ENST00000260682.7_hg38_4_9 54 1 0 chr10:94947779-94947939+ [161-214]
ASPCDPTFILGCAPCNVICSIIFHKRFDYKDQQFLNLMEKLNENIKILSSPWIQ
>ENST00000260682.7_hg38_5_9 59 0 0 chr10:94949108-94949284+ [215-273]
ICNNFSPIIDYFPGTHNKLLKNVAFMKSYILEKVKEHQESMDMNNPQDFIDCFLMKMEK
>ENST00000260682.7_hg38_6_9 47 0 1 chr10:94972104-94972245+ [274-320]
EKHNQPSEFTIESLENTAVDLFGAGTETTSTTLRYALLLLLKHPEVT
>ENST00000260682.7_hg38_7_9 63 1 0 chr10:94981183-94981370+ [321-383]
AKVQEEIERVIGRNRSPCMQDRSHMPYTDAVVHEVQRYIDLLPTSLPHAVTCDIKFRNYLIPK
>ENST00000260682.7_hg38_8_9 47 0 1 chr10:94986033-94986174+ [384-430]
GTTILISLTSVLHDNKEFPNPEMFDPHHFLDEGGNFKKSKYFMPFSA
>ENST00000260682.7_hg38_9_9 61 1 0 chr10:94988847-94989028+ [431-491]
GKRICVGEALAGMELFLFLTSILQNFNLKSLVDPKNLDTTPVVNGFASVPPFYQLCFIPVZ

And these nucleic acid sequences:
>ENST00000260682.7_hg38_1_9 168 0 0 chr10:94938683-94938850+
ATGGATTCTCTTGTGGTCCTTGTGCTCTGTCTCTCATGTTTGCTTCTCCTTTCACTCTGGAGACAGAGCTCTGGGAGAGGAAAACTCCCTCCTGGCCCCACTCCTCTCCCAGTGATTGGAAATATCCTACAGATAGGTATTAAGGACATCAGCAAATCCTTAACCAAT
>ENST00000260682.7_hg38_2_9 163 0 1 chr10:94941858-94942020+
CTCTCAAAGGTCTATGGCCCTGTGTTCACTCTGTATTTTGGCCTGAAACCCATAGTGGTGCTGCATGGATATGAAGCAGTGAAGGAAGCCCTGATTGATCTTGGAGAGGAGTTTTCTGGAAGAGGCATTTTCCCACTGGCTGAAAGAGCTAACAGAGGATTTG
>ENST00000260682.7_hg38_3_9 150 1 1 chr10:94942192-94942341+
GAATTGTTTTCAGCAATGGAAAGAAATGGAAGGAGATCCGGCGTTTCTCCCTCATGACGCTGCGGAATTTTGGGATGGGGAAGAGGAGCATTGAGGACCGTGTTCAAGAGGAAGCCCGCTGCCTTGTGGAGGAGTTGAGAAAAACCAAGG
>ENST00000260682.7_hg38_4_9 161 1 0 chr10:94947779-94947939+
CCTCACCCTGTGATCCCACTTTCATCCTGGGCTGTGCTCCCTGCAATGTGATCTGCTCCATTATTTTCCATAAACGTTTTGATTATAAAGATCAGCAATTTCTTAACTTAATGGAAAAGTTGAATGAAAACATCAAGATTTTGAGCAGCCCCTGGATCCAG
>ENST00000260682.7_hg38_5_9 177 0 0 chr10:94949108-94949284+
ATCTGCAATAATTTTTCTCCTATCATTGATTACTTCCCGGGAACTCACAACAAATTACTTAAAAACGTTGCTTTTATGAAAAGTTATATTTTGGAAAAAGTAAAAGAACACCAAGAATCAATGGACATGAACAACCCTCAGGACTTTATTGATTGCTTCCTGATGAAAATGGAGAAG
>ENST00000260682.7_hg38_6_9 142 0 1 chr10:94972104-94972245+
GAAAAGCACAACCAACCATCTGAATTTACTATTGAAAGCTTGGAAAACACTGCAGTTGACTTGTTTGGAGCTGGGACAGAGACGACAAGCACAACCCTGAGATATGCTCTCCTTCTCCTGCTGAAGCACCCAGAGGTCACAG
>ENST00000260682.7_hg38_7_9 188 1 0 chr10:94981183-94981370+
CTAAAGTCCAGGAAGAGATTGAACGTGTGATTGGCAGAAACCGGAGCCCCTGCATGCAAGACAGGAGCCACATGCCCTACACAGATGCTGTGGTGCACGAGGTCCAGAGATACATTGACCTTCTCCCCACCAGCCTGCCCCATGCAGTGACCTGTGACATTAAATTCAGAAACTATCTCATTCCCAAG
>ENST00000260682.7_hg38_8_9 142 0 1 chr10:94986033-94986174+
GGCACAACCATATTAATTTCCCTGACTTCTGTGCTACATGACAACAAAGAATTTCCCAACCCAGAGATGTTTGACCCTCATCACTTTCTGGATGAAGGTGGCAATTTTAAGAAAAGTAAATACTTCATGCCTTTCTCAGCAG
>ENST00000260682.7_hg38_9_9 182 1 0 chr10:94988847-94989028+
GAAAACGGATTTGTGTGGGAGAAGCCCTGGCCGGCATGGAGCTGTTTTTATTCCTGACCTCCATTTTACAGAACTTTAACCTGAAATCTCTGGTTGACCCAAAGAACCTTGACACCACTCCAGTTGTCAATGGATTTGCCTCTGTGCCGCCCTTCTACCAGCTGTGCTTCATTCCTGTCTGA

Associating the pocket protein amino acid numbers in exons 2,3,5,7,9 with genomic sequence using these FASTA sequences is left as an exercise for the reader.

source: Stack Exchange

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

sort: votes


text:
filter out FASTQ reads which are shorter
I tried to filter out FASTQ reads which are shorter than 259 bp with
bioawk -cfastx 'length() >= 259 {print "@" " " "\n""\n+\n"}' good/SZ005_NoIndex_L002_R1_009.good.fq.gz \
  | gzip >  good-filtered/SZ005_NoIndex_L002_R1_009.good-filtered.fq.gz`.

Unfortunately, the output was wrong:
@

+

@

+

@

What did I miss?
Thank you in advance,

If you have paired-end reads, this solution keeps the two files in-sync (i.e. discard pairs where one of the two reads is shorter than 259). Also, it uses only Unix tools without requiring external programs:
paste <(zcat S1_L001_R1_001.fastq.gz | paste - - - -) <(zcat S1_L001_R2_001.fastq.gz | paste - - - -) \
| awk -v FS='\t' 'length($2) >= 259 && length($6) >= 259' \
| tee >(cut -f 1-4 | tr '\t' '\n' | gzip > r1.fq.gz) \
| cut -f 5-8 | tr '\t' '\n' | gzip > r2.fq.gz

You need to use variable $seq. For example:
bioawk -cfastx 'length($seq)>=259{print "@"$name"\n"$seq"\n+\n"$qual}' test.fq.gz

There are also more convenient tools like seqtk and seqkit. With those, you may
seqtk seq -L260 test.fq.gz

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


text:
Finding transposable elements using RepeatMasker
I'm using RepeatMasker to detect, classify the Transposable elements. My Input is a eukaryotic non-reference genome.
I made a run via RepeatMasker many times to Mask the TEs, but return 0 Annotation tables.
Further, I used a different -species option each time, but Annotation results contain only short and simple repeats.
Do I have to create a unique library to detect them?
The command I used:
./RepeatMasker -no_is -noint - species mammal Genome.fna

Results:
sequences:         32573
total length: 2004063690 bp  (1981588036 bp excl N/X-runs)
GC level:         41.28 %
bases masked:   33282890 bp ( 1.66 %)
==================================================
               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:         0            0 bp    0.00 %

Total interspersed repeats:        0 bp    0.00 %

Small RNA:            0            0 bp    0.00 %

Satellites:           1          267 bp    0.00 %
Simple repeats:  663679     27953718 bp    1.39 %
Low complexity:  123705      6154714 bp    0.31 %
==================================================

is there any alternative or suggestions for a tool that accomplishes this task, scan the sequence and detect these Elements

Welcome to Bioinformatics. Could [edit] the question to add your comment to the body of the question?  Also could you clarify what do you mean with different species and options ? Are you studying a mammal organism or not?

Maybe use something like RepeatScout to build a library for RepeatMasker, relevant https://bioinformatics.stackexchange.com/a/347/104

looks like this can a solution for your problem
https://www.biostars.org/p/154290/

what I meant by different species . its an option in Repeatmasker to look under a specific range of species. I'm studying eukaryotic sequence so I thought mammal species could be the closest.

does the Transposable elements "equals" the repeats that give by the RepeatMasker, which also can be download by the UCSC table browser.

That's one answer ^, not sure if you NEED to do that, which wouldn't actually classify your elements as LTRs, it would just find repetitive sequences in your new genome (as you discovered).
To get them classified one standard approach is to use evidence of homology to sequences of LTR elements with classifications already, and you can get those in Repbase. If you want to download Repbase, you have to sign up for an account with GIRI but I think they give to anyone who asks. Then you can download the section of Repbase you're interested in (e.g. the whole thing, LTR, DNA-TE, ...) and provide the path to the downloaded library (in FASTA format I think) using -lib, like
RepeatMasker -lib Repbase_LTR.fa genome.fa
Then your RepeatMasker output table should have info about other kinds of TEs.
If you want to check out my approach to classification by homology you might find some ideas you'd like to try. I used Dfam+nhmmer and Repbase+tblastx to identify evidence for LTR retrotransposon classification.
My software pipeline is available as open-source software. If you want detailed high-quality annotations of LTR retrotransponsons, PhyLTR will automatically annotate putative protein-coding sequence regions, classify them, and remove false-positives like a tandem array of DNA TEs which look like LTR-Rs (good diagram of LTR-R false positives here: https://github.com/oushujun/LTR_retriever/blob/master/Manual.pdf)
Phylogenetic Analysis of LTR retrotransposons
https://github.com/mcsimenc/PhyLTR

This is indeed an interesting approach. Do you have any comments on the actual homology or sequence conservation of different repetitive element "entities" across multiple species? Also wouldn't the approach you mention miss non-LTR TEs? LINE-1 elements are the most abundant TEs in humans but not have LTRs.

I'm not sure what you're thinking of exactly, but I can say that in the all plant genomes I have annotated using this approach resulted in clusters of "Gypsy" "Copia" elements where domain content and organization independently affirmed their classification as such. Also some lower-level "taxa" like Ogre elements were recognizable by domain architecture. This work is meant for LTR elements only. It would work fine for metazoan endogenous retroviruses and LTR retroviruses, but not for LINEs. I think the author of the post mentioned LTRs in the original post on stackoverflow.

Dears : i appreciate these knowledge, alight me towards my projects , i have gone far in my analysis using ltrharvest + ltrdigest using 314 protein profiles from Gydb database and trna library obtained using tRnascan, then filtration candidates with protein hits only Remain, now i do have the coordinates and the sequences of the TE class1 , i m trying to classify them into subfamilies ex (copia,gypsy...etc) but i failed , as well as class 2 DNA, MITES and Helitron Results obtained with different pipeline but fail to classify them into subfamilies @matt

@story my project is to report all the Transposable Element in my Draft Genome Class1 LTR and NON-LTR . Class 2 Dna transposons , successfully the LTR structure easily can be detected by Denovo tools , but non-LTR must of the tools need windows operating system like MegaScan non-ltr any suggestions

@matt i did registered to download the database but couldnot find the way to specify only LTR or Dna it comes all together . i tried /// tBlastn - query LTR-filtered_digest.fasta - subject mem.ref -out results ///  vrtb.ref its one of many files downloaded from repbase refer to vertebrate but no hits are found

@emaz, sounds good, nice work. The GIRI Browse page, https://www.girinst.org/repbase/update/browse.php, is place I was thinking for downloading the kinds of repeats separately. Other than that you might need to do some manual steps or scripting to mold those files how you want, e.g. cat mem.ref vrtb.ref >> repbase.fasta to join two text files. If I had the compute power I might look at all of them together. If I remember correctly, RepeatMasker gets taxon info/descriptors from the FASTA sequence header/ID by default. MGEscan-nonLTR - https://mgescan.readthedocs.io/en/latest/nonltr.html.

@emaz, tblastx is the one to use, I believe. tblastn searches a translated nucleotide database (=protein seqs) using a protein query -- If you give nucleotide FASTA as query input it will still treat it like a protein, i.e. ACG = Alanine-Cysteine-Glycine. tblastx will ultimately do a protein protein search just like tblastn but it treats the query input as a DNA sequence and translates it for you. I'm curious if you will get some hits using tblastx. tlblastx is one of the slowest BLASTs because it translates both the query and subject (database) sequences in all six reading frames.

^so tblastx is one of the most sensitive BLASTs to the kinds of haphazard sequence mutations that tend to happen in loci undergoing neutral evolution, like probably many TEs are, mostly; i.e. there will probably be a lot of frameshifts and the nucleotide seqs will change much faster than their encoded protein seqs. tblastx can find fragments of similar encoding sequences, which makes it helpful for comparing fragmented TE sequences, which most are.

im so glad to have you here @matt so as i understand will run // tblastx -query ltr-sequences.fasta -subject mem_vert_combined.ref //  to classify them into subfamilies ?

great, I'm glad to help. What your ^ command will do is run tblastx, which if you add -out OUTPUTFILENAME and -outfmt 7, will create a potentially very large tab-delimited file with information about which parts of each sequence have matches (high-scoring sequence alignments with) your input database sequences. You will probably need to do some manual scripting at this point to decide and add your classifications to your TE file(s). I have some Python scripts I have made for doing this kind of thing; if you get there and want some more help get a hold of me and I can send you the scripts.

Let us continue this discussion in chat.

You need to model the repeats in your de novo genome.
See:
http://www.repeatmasker.org/RepeatModeler/ and
https://www.biostars.org/p/154290/ (from previous answer)
You first build a database named "name_of_your_database" (insert whatever you want to call it) and run it on the fasta file of your genome (e.g. your_genome.fasta).
Make you find the right paths for the software.
## build the database
BuildDatabase -engine ncbi -name name_of_your_database your_genome.fasta

## run the modeler (in this case using ncbi)
RepeatModeler -database name_of_your_database -engine ncbi

In case someone finds this years later: OP's problem was the use of the "-noint" flag. It stands for "no interspersed", so the OP explicitly asked RepeatMasker  NOT to look for transposable elements. If they omitted the flag, RepeatMasker would have automatically downloaded mammal TE sequences (nowadays it uses Dfam). Even if they had a custom library like suggested above, they would still need to delete "-noint". Check your commands!


text:
How to combine multiple kernels of large sample datasets?
I have multiple large sample datasets in matrix format (each has 15000 rows and 5-50 columns) corresponding to different experiments. Each matrix contains the same number of samples(rows) but the variables(columns) are not the same. My objective is to cluster the samples on the basis of all the experiments.
I tried to use Unsupervised multiple kernel learning (UMKL) to integrate the datasets followed by kernelPCA using "mixkernel" package in R (https://www.ncbi.nlm.nih.gov/pubmed/29077792). The UMKL step calculates kernels for each dataset and combines the kernels using 3 different approaches: 1) calculating a consensus kernel from multiple kernels 2) calculating a sparse kernel preserving the original topology of the data 3) calculating a full kernel preserving the original topology of the data
The kernel calculation step was fine but the kernel integration step (all 3 approaches) runs very long and my computer hangs.
Is there any way to handle this problem? More specifically, is there any way to handle multiple kernel integration for large sample datasets?
Any suggestion alternative to using kernel methods will also work.

Welcome to Bioinformatics Stackexchange @SD1024. A recently published algorithm known as MOFA (Multi-Omics Factor Analysis, paper, github) is generating a lot of interest, and is designed answer exactly the sort of biological question you are describing. It claims to extract axes from multiple matrices with overlapping samples, but not necessarily overlapping features and handle missing data elegantly. I havn't tried it myself you, but I'm quite excited to do so.

Thanks for the suggestion. Actually, I used that package. If I have understood correctly, MOFA extracts the latent variables from individual data modality nicely but such variables do not integrate multiple modalities together. For example, If multiple Chip-Seq and CpG methylation assays are done on multiple genomic sites and across multiple patient samples/cell lines, MOFA extracts one set of latent factors for ChIP-seq and one set of latent factors for CpG methylation. It does not extract any latent factor (or a set of latent factors) that represents both the assays simultaneously.

As I say, i've not actaully used it myself, but the paper definately claims that the factors intergate multiple assays. For example, see their figure 2B, where the plot shows contributions to factor 1 from all assay types, but only contributions from drug response, mutation data and transcription (i.e. not from methylation) to factor 2

You are right. It is described in equation 1 of materials and method. A great help. Thanks a lot. :)

You can use regularized general canonical correlation analysis to cluster those matrices. I would suggest RGCCA, which allows to do this for an arbitrary number of matrices but it  doesn't use kernels.
There has been some efforts to use kernels on the canonical correlation analysis but there isn't any implementation (yet), you can read more on the related literature, specially the Kernel Generalized Canonical Correlation Analysis.
However in this kind of analysis you will need to come up with a design of how those matrices are related between them. It is not like a PCA where you give some data and you get some data back that summarize the matrix.


text:
python package for NNI neighbors
I am working on protein sequence data files for reconstructing phylogenetic tree and I need to generate all NNI-neighbours of the tree (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation). Is there any python package/sorce code I can access to, which can retrieve all NNI neighbors of a tree.
Earlier I have used biopython 1.72 to generate tree from my dataset. But for now I'm looking for a python package which can do NNI search on my tree and return all NNI-neighbors.

Both the ETE Toolkit and DendroPy offer methods of node removal and reattachment - I would think it would be easy enough to hack together an NNI movement function with them.
You could also investigate the source code of BioPython's NNITreeSearcher() function.

I have used biopython among them, but my concern is to return all the trees which are reconstructed before the best tree is reached, and biopython NNITreeSearcher() results one best scored tree.

That is why I suggested searching the source code to see if you can find the code controlling the NNI moves and use that.

I have studies biopython module and I found that _get_neighbors() function generates all NNI neighbor trees of the given tree.


text:
Visualize network (graph) using "line" layout
I want to visualise a network in a way that the nodes are along a line (and the edges are not on that line).
I have included a simple example below.

Which network visualisation programs (Cytoscape, Gephi, R igraph etc) do have this layout option?

That plot is called an Arc diagram
you can google that and use the package or program that best suits you.
Here are 2 examples of the use of Arc diagrams in R:
Example 1:
Example 2:

I'm the author of gravis, an open source package for interactive network visualization in Python. It supports many ways to influence visual properties with data, such as node size, color or location. It also allows to control edge curvature. Therefore it's possible to place all nodes along a line (or another form) and have edges with a chosen degree of curvature between them.

Example 1

Example 2

Edit: Python code that reproduces your image with NetworkX and gravis
import gravis as gv
import networkx as nx

g = nx.Graph()
g.add_node(1, x=0, y=0)
g.add_node(2, x=50, y=0)
g.add_node(3, x=100, y=0)
g.add_node(4, x=150, y=0)
g.add_edge(1, 2)
g.add_edge(1, 3)
g.add_edge(3, 4)
g.graph['node_color'] = 'blue'
g.graph['node_size'] = 15
gv.d3(g, edge_curvature=1.5, show_node_label=False)

Result inside a Jupyter notebook (optional, just for convenience):


text:
Aligning sequence and comparing it against primer
I am looking to show how a primer is consistent among some genomic data. I have a primer of about 23bp and looking to compare it to about 5000 genomic sequences of 10kb. I am unsure how to format it the way I need to, since I am just trying to show the nucleotide differences of the primer to the genetic sequences. Here is what I need:
>     1. Cut out the area that my primer is located and about 20bp down each end.
>     2. Show only the bases that are different from my primer in my analysis.
>      ex: Primer: -----------ATGTGGAAGCAAATATCAAATGA---------
>          Gene:   ATGACCATACG----C--------------T---ATCGTAGGG
>                  ATGAGCATACC-----A----T--------T---TTCGAACGC

The data I am using is all dengue sequences (all serotypes) and the primer with the following code: ATGTGGAAGCAAATATCAAATGA.
I am using Biostrings, msa, and seqinr. I use ncbi to get the genetic sequences and using FASTA files.
Thanks!
Edit: Solution:
I used MAFFT to align the data, working very efficiently and put it into a fasta file. Then I cut the nucleotides I did not need by:
1. Finding the part i wanted and cutting 10 nucleotides before it and after it.
2. Putting only the [start:end] into a new list, which I will then manipulate.
However, I am still having difficulty finding a way of showing my data nicely because I can't find a suitable package. Any advice?
Edit:
Hi all,
Just wanted to show the code I made to solve the problem.
    #So the we will be using a dataframe that has the sequences as a string
    #Then I created a temporary dataframe to store each sequence, with the primer at the top

    tempDF <- data.frame(matrix(primer1, nrow = 1, ncol = 1),
    stringsAsFactors = FALSE)

    for(i in sequenceAsString){
      #Temporary vector for each sequence after compared to primer.

      temporary1 <- ""
      temporary1 <- paste(temporary1, substring(i,1,11), sep = "")

      for(j in 12:(nchar(i)-20)){
        if(substring(i,j,j) == substring(toString(reversePrimer1), j-12,j-12)){
          temporary1 <- paste(temporary1, "-", sep = "")
        }else{
          temporary1 <- paste(temporary1, substring(i,j,j), sep = "")
        }
      }
      temporary1 <- paste(temporary1, substring(i,nchar(i)-19,nchar(i)), sep = "")
      tempDF <- rbind(tempDF,temporary1)
    }

So, I essentially take the sequence I chose apart and compare each character to that of the reverseprimer. I then replace the ones that are the same with an '-'
Thanks everyone for their help!

“Since that is too much for my computer to do” — No. Performing a short pairwise alignment like this is trivial. Running it against 5000 sequences of length 10knt should take a few seconds at most. Maybe a bit longer in R using pairwiseAlignment but totally feasible. And if performance is really of the essence, you can implement e.g. Myers Bit Vector since the primer is so short, and the whole search could happen in sub-second time. Either way, don’t perform a multiple sequence alignment.

The issue you are having is not entirely clear. When you say "However, it was difficult because to accurately predict if you would need to have it aligned.": Does this mean you do not know how to select the homologous region across all genes that encompass the primer? You may want to be more specific, ie provide detailed example, what you tried and exactly where issue is.

Thanks a lot Konrad and Vince, very helpful. I am clearly just learning the ropes! I will get back with the result.

So, I have figured it out using a combination of MAFFT and R. I still have not found a way to present the data in the format I mentioned, might you have any examples of packages that can be used?

If I understand you correctly, after your most recent edit, you have 5,000 sequences that are 63 bp long and in each you simply want to replace the nucleotides from 21 through 43 with - unless that nucleotide does not match the primer. Since - is normally used to represent gaps, I recommend using . instead for matching nucleotides.
I know you're looking for an R solution, but I don't usually use R for string manipulations and don't have time at the moment to figure it out. So, assuming this might be urgent for you, below is a Python solution. This uses your example data which shows 11 nucleotides before the primer and 9 after it (so you'll need to adjust it to work with your actual case in which 20 nucleotides are before and after the primer). In Python, strings are arrays of characters, by the way.
Here are your example sequences are in the file seqs.txt:
ATGACCATACGATGTCGAAGCAAATATCAATTGAATCGTAGGG
ATGAGCATACCATGTGAAAGCTAATATCAATTGATTCGAACGC

Here is the code that gives you the expected output:
#!/bin/env python
import sys
primer = 'ATGTGGAAGCAAATATCAAATGA'

print('-' * 11 + primer + '-' * 9 )
with open('seqs.txt','r') as f:
    for gene in f:
        sys.stdout.write(gene[0:11])
        for gene_index in range(11,34):
            primer_index = gene_index - 11
            if gene[gene_index] == primer[primer_index]:
                sys.stdout.write('.')
            else:
                sys.stdout.write(gene[gene_index])
        print(gene[34:])

Running this gives you:
-----------ATGTGGAAGCAAATATCAAATGA---------
ATGACCATACG....C..............T...ATCGTAGGG
ATGAGCATACC.....A....T........T...TTCGAACGC

Thanks a lot Christopher, I really appreciate that you took the time to help me out. I will definitely use this and potentially try to implement it in R. Thanks!


text:
Sequence alignment using Markov Model
I am learning about applying Markov model to sequence alignment. The prof says that the transition probabilities from a gap-residue alignment to a residue-gap alignment and vice versa are both 0. Is there any biological/mathematical reason behind this statement? Why are the (X,Y) and (Y,X) cell 0? This is a lecture slide of lecture 1, week 4 of the "Bioinformatics: Introdcution and Methods" course on coursera.

If I understand your question correctly, then I think for case of pairwise alignment, there is a simple explanation.
I believe the key insight is that: a mismatch should always score better than a gap.*
This follows biologically since the insertion/deletion (indel) rate is roughly 1/10th that of the substitution rate (i.e. the occurrence of single nucleotide changes), at least in vertebrates. (This varies across the tree of life but I think the substitution rate virtually always exceeds the indel rate.)
To understand why this matters, consider an example:
ATG-AG
ATGT-G

This is an 'impossible alignment' under the probabilities you gave since here we have a transition from a gap-residue alignment to a residue-gap.
However, under our assumption that mismatches are more likely biologically than indels, the correct alignment should be:
ATGAG
ATGTG

Indeed, the latter does look like a better alignment.
This also follows for more complex examples, so this:
ATG--AAG
ATGTT-AG

Becomes this:
ATG-AAG
ATGTTAG

(Or this:
ATGA-AG
ATGTTAG

)
* Strictly, I mean a substitution should score better than an indel (with the associated gap opening and extension penalties). In fact, for the assumption to always be true, a run of mismatches should still score worse than a single indel. This may not always be a correct assumption, consider this example below, is the true alignment case 1) or 2) or something else? Or is in fact a global alignment bad here and this should be split into 2 local alignments? Is there a likely biological mutational event that could explain this? I ask these questions just to point out it is not black-and-white, I don't have clear answers
1)
CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

2)
CGTACGTAGAGGAATGCCCCCCCCCAGCAACGTAGCAT
CGTACGTAGAGGAATGTTTTTTTT-AGCAACGTAGCAT

All Chris_Rands said is correct: you set the probability of $X\to Y$ and $Y\to X$ to 0 to forbid adjacent insertions/deletions in the alignment. A lot of textbooks including some classical ones use this rule, but in fact, the rule is questionable. It is easier to see this from Smith-Waterman alignment under the affine gap penalty, which is largely the non-probabilistic view of paired HMM.
With the affine gap penalty, a gap of length $k>0$ is scored as
$$
g(k)=-(d+k\cdot e)
$$
where $d\ge0$ is the gap open penalty and $e>0$ is the gap extension penalty. Suppose we are using a simple scoring matrix where a mismatch gets $-b$, $b>0$. We may see an insertion immediately followed by a deletion (and vice versa) if $b>2e$. It is actually not so difficult for this to happen. For example, for the human-mouse alignment (see the blastz paper), $e=30$ and $b$ is ranged from 31 to 125. It is possible that an $X\to Y$ transition is preferred in the alignment.
Theoretically speaking, it makes more sense to consider immediate transitions between insertions and deletions. In practice, though, the difference between allowing/disallowing such transitions is probably minor most of time.
EDIT: on Chris' example
CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

If we use a scoring matrix with $b>2e$ but disallow adjacent ins-to-del transitions, we will probably end up with an alignment like
CGTACGTAGAGGAATGCCCCCCCCC-------AGCAACGTAGCAT
CGTACGTAGAGGAATG--------TTTTTTTTAGCAACGTAGCAT

This alignment score will be lower.


text:
Curl error when installing Bioconductor packages with Microsoft R Open 3.5.0
I'm trying to install Bioconductor packages with Microsoft R Open 3.5.0 seems to fail not sure what is the issue.
This is what I'm trying to do
source("https://bioconductor.org/biocLite.R")
biocLite("BiocInstaller")

But I keep getting an error about "curl" had non-zero exit.
How can I install the packages of Bioconductor?

Please [edit] your question and clarify. Your first paragraph seems to be about installing R but the last one suggests you have already installed it. Please don't ask multiple questions in a single post, [edit] and ask one, specific and clear question. Thanks!

updated my question ..iim getting this "curl' non-zero exit which wanted to post the error but closed my terminal so lost the exact error message

Then please run the command again and then add the error message to the your question. We can't help you solve an error if you don't show it.

I can't comment yet but from what you said it seems that you need to install additional libraries like the one you mentioned: "curl".
from the R terminal you should type:
install.packages("curl")

then R will install the package. Not sure if there are some other dependencies. If yes R will return you some error message in which there are some package names you need to install.
EDIT:
of course, after installing all the required packages you should be able to install biocLite too.

cd /opt/microsoft/ropen/3.5.1/lib64/R/library/  inside the path i can see the curl so why still i get the error i'm not sure

try loading the curl package from the R terminal and see if it actually works (library(curl)). Maybe some libraries are missing. I use linux and curl requires additional libraries that are not installed by R but are system's libraries.

yeah i did check the system level libraries as i can see using this "locate libcurl "so still not sure why the non-zero exit

finally after some time it did work.

sorry for not replying! if you solved the issue just post how you did it. This will help other people having the same problems!

earlier it was not showing any error no idea what was going on but later it showed up the dependencies  error which i installed it sort of fixed the problem ..but the way i used to update using CRAN is now no longer working with Microsoft R open here i can;t define the repository how to fix that any suggestion ?

@krushnachChandra, sorry for this very late reply, I must have lost your comment...have you solved your issue? otherwise I believe you should open another question with more details.

well thank you for your late response , i had few other issues as well .So i reverted to the base R


text:
What is the correct way of dealing with the analysis of data from flow cytometry?
I would like to detect a change in expression of a molecule present on a cell type by flow cytometry.
Assuming I am able to detect, using an antibody, a signal that represents the amount of the molecule I'm interested into. Also assume that the same experiment will have some noise (ideally near to zero) represented by an isotype antibody. What if the noise is somehow higher than a minimal threshold (i.e. = 1%)? How can I deal with the noise? How can I compare different samples if the noise is higher than the threshold?
I was thinking about removing those samples that have a noise above the 1% but is this ok, from a statistical point of view?
Another possibility is: can I subtract the noise from the signal and use this information for the final statistical analysis? Is this correct? Is it possible to subtract a signal from another? How is this done?
EDIT just to clarify a bit:
I have a pool of neutrophils which, in resting condition, do not express the molecule on their surface. When stimulated, these cells are supposed to show an increased amount of this molecule (again, on their surface). We would like to quantify the change in the amount of molecule (we used FACS).
We are able to quantify the quantity of surface-molecule but we also get some aspecific signal (the noise due to the isotype antiboby). What we try to do is to keep the isotype signal lower than 1% using a gate-based methodology. What if the isotype signal we get is higher than the threshold (i.e. >1%)? The antibody sometimes work some others don't. What should we do when it returns too much noise? Should we discard the data or should we perform a normalisation over the surface-molecule using the noise?

If your goal is to simply "to quantify the change in the amount of molecule," then you can do that at the population-average level by simply taking the average antibody signal for the different cell populations (stimulated vs. not stimulated). In that case, you don't have to throw out samples with noise above a threshold - since they are still part of the population and so should be used when calculating average.
If you want to quantify "the change in the amount of molecule" at the single-cell level, then things are a bit more complicated, and the best answer depends in part on things like the shapes of the distribution of signal from the stimulated and non-stimulated populations.
And related to that, when you say "noise above the 1%," do you mean that the Coefficient of Variation (CV) is greater than 1%? Or something else?


text:
Error in annotating row names as data points in a PCA plot
I made a PCA plot and was trying to plot the eigenvenctors in R so that each data point is actually the sample name. All the samples are in column 1. The following are the R commands I am using:
    fn <- "Plot.txt"
    evecdat <- read.table(fn, col.names=c("Sample", "PC1", "PC2", "PC3", "PC4", "PC5", "PC6", "PC7", "PC8", "PC9", "PC10", "pc11","pc12", "pc13", "pc14", "pc15","pc16", "pc17", "pc18", "pc19", "pc20", "pop", "color"))
    jpeg("Plot1.jpeg", width = 20, height = 20, units = 'in', res = 300)
    identify(evecdat$PC1, evecdat$PC2  = NULL, labels = seq_along(evecdat$Sample), pos = TRUE, n = length(evecdat$Sample), plot = FALSE, atpen = FALSE, offset = 0.5, tolerance = 0.25, …)

    plot(evecdat$PC1, evecdat$PC2, col=as.vector(evecdat$color), xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex=1.5, cex.axis=1)
    legend("topright", ncol=1, legend=c("X", "Y", "Z"), col=c("Blue", "Red", "Green"), pch="24", cex=2, bty = "o")
    while (!is.null(dev.list())) dev.off()

Here, all the row names are under the Sample column name that I defined earlier.

We don't have the "Plot.txt" data file so we can't test, but what problem do you have with the code? You get the plot without the name of the samples? Have you looked at  text?

Maybe try
plot(evecdat$PC1, evecdat$PC2, xlab="PC1", ylab="PC2", pch="24", cex.lab=2.5, cex.axis=1, type="n")
text(evecdat$PC1, evecdat$PC2, evecdat$Sample, col=as.vector(evecdat$color))

This should plot the sample names as text (taken from the evecdat$Sample column), colored as specified in the evecdat$color column. See ?text for details.


text:
dog coordinates (canFam3) to human coordinates (hg19)
I've converted dog coordinates to human using UCSC LiftOver. These are 200bp intergenic regions that are differentially methylated from normal dogs to cancer dogs. I've converted these to human coordinates and found that a lot of them overlap with differentially methylated regions we found in the human model.

Is this okay to do?
What stringencies should I check or modify, the min ratio of bases needed to overlap?
Lastly, is LiftOver taking the 200bp dog sequence and saying that at least 10% of the sequence has to align with a region within the human genome?

I could not find an answer online.

Based on this old answer, this is not recommended for liftOver. For between-species liftovers, it is claimed that you are likely to want a different tool such as pslMap. For a lot of rather detailed discussion of how liftOver operates in this context and how to parameterize properly, see here.
See also this old Biostars question addressing the issue of mapping not between versions of the same assemblies but different genomes entirely. This thread also notes that using liftOver can be dangerous. That isn't to say that you can't learn something interesting, just that you might make mistakes.
pslMap apparently drops non-syntenic regions, which may lose a lot of data, but it is supposedly making rather fewer assumptions about assembly relationships.
orthoMap is apparently a similar tool to pslMap, but I can't find a way to run it except for a very old and restricted web interface. There are other tools with very similar names ("orthomap") but somewhat different uses, which may increase confusion.
A counterpoint to all the foregoing is that this rather detailed sum-up of liftover tools does not seem to be very concerned about the issues of using tools like liftOver between species. It specifically uses dog to human liftover as the example, just as you want to. It's notable that the assumption is made in that page that a high-quality chain file exists, which does seem to be true for you.
To go through your questions:

depends, see above.
For suggested parameters for running blat/chain generation/liftOver for this application, see here.
This is derived in part from the [chain file] I believe rather than liftOver parameters (http://genomewiki.ucsc.edu/index.php/LiftOver_Howto). However overall this is set by the parameters you used in generating a chain file, which probably boils down to parameters of the workflow here. liftOver itself is just relying on the chain alignments for each region with the best score.


text:
Reproducing population structure from 1000 Genomes Data
I recently tried to reproduce population structure with 1000 Genomes Project data (Phase 3, 2504 people).
According to supplementary material,

To reduce the dataset, we used VCFtools to only keep biallelic, non-singleton SNV sites that were a minimum of 2 KB apart from each other.
  These variants were then merged and converted into PLINK format, giving us 1,286,213 throughout all the autosomes.
  PLINK was then used to filter the variants further, keeping only variants with a minor allele frequency of at least 0.05,
  and to convert the variants into binary PED format. After this filter 193,634 sites remained with data from 2,504 individuals in 26 populations.

And from ftp repository, I found the exact command-line to do it.

To get biallelic, non-singleton SNVs, all minimum of 2KB apart from each other

./vcftools --gzvcf
ALL.chr1.phase3_shapeit2_mvncall_integrated_v4.20130502.genotypes.vcf.gz
--thin 2000 --min-alleles 2 --max-alleles 2 --non-ref-ac 2 --plink --chr 1 --out ALL.chr1.phase3_shapeit2_filtered.20141217

To get variants with minor allele frequency of at least 0.05

./plink --file ALL.wgs.phase3_shapeit2_filtered.20141217.ped --out
ALL.wgs.phase3_shapeit2_filtered.20141217.maf0.05 --make-bed --maf 0.05

To get population structure with ADMIXTURE program (v1.23)

./admixture ALL.wgs.phase3_shapeit2_filtered.20141217.maf0.05.bed K

So I applied the same command (although versions of programs were not exactly the same) to the phase 3 data, and got different results (top: from supplementary data, bottom: my reproduction)

I basically did the same thing as in the supplementary document:

(for each chromosome \$ch) bcftools view -v snps -O v $gzvcf | vcftools --vcf - --thin 2000 --min-alleles 2 --max-alleles 2 --non-ref-ac 2 --plink --chr ${ch} --out 1KGP_v5a_20130502.chr${ch} [bcftools-1.9, vcftools 0.1.15]
(for each chromosome \$ch) $plink --file 1KGP_v5a_20130502.chr${ch} --maf 0.05 --make-bed --out 1KGP_v5a_20130502.chr${ch}.maf0.05 [plink 1.9b]
(merge)
admixture 1KGP_v5a_20130502.ALL.maf0.05.bed 5 [ADMIXTURE 1.3]

I would like to know if there's some step I've missed on the process. Otherwise, I cannot understand why the population structure look so different.

Welcome to the site. As you mention that you use different software version: What versions are you using? Could you try the same version to see if you are able to reproduce the same figure?

Did you perform your steps on all the chromosome or just chromosome 1?  Also in step 1 why did you skip the bedtools command that was used in the original paper?

@Llopis : I've noted the software versions I used at the end of command lines. When started with the same .bed/.bim files from the [repository] (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/admixture_files/), I was able to get the same output with ADMIXTURE 1.3 as with ADMIXTURE 1.23 (original paper). So, at least, ADMIXTURE version doesn't seem to matter.

@Bioathlete : I did for all autosomes. I followed the steps in the [file] (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/admixture_files/README.admixture_20141217) & supplementary section 9.7. I thought those steps can be directly applied to the released VCFs. Can you point me the bedtools command you mentioned?

sorry I meant bcftools not bedtools.  bcftools view -v snps -O v $gzvcf

@Bioathlete : The commands at the bottom were mine. Those at the top of the figure are from the supplementary document. Sorry about the confusion.

Ok then the question is reversed but still valid.  1) Why did you use bcftools? 2) Is the image from the paper using the whole genome or just chromosome 1?

@Bioathlete : 1) I used bcftools to select SNPs only, because it is described as they used SNVs. But the 1000G phase 3 VCF contained indels & SVs. 2) Image in the paper was also whole genome. After all, paper used 193,634 sites (over whole autosome) as input to ADMIXTURE, while I used 199,657 sites.

Did you check how many iterations they ran ADMIXTURE for, and how many you did? It may have yet to converge.


text:
AttributeError: module 'scater' has no attribute 'normalize'
I am using rpy2 to run normalization using scran package. After calling computeSumFactors I want to do the actual normalization calling normalize() function from scater package:
def normalize(counts):
utils = importr('utils')
utils.chooseCRANmirror(ind=1)

packnames = ('Matrix', 'limSolve')
names_to_install = [x for x in packnames if not isinstalled(x)]
if len(names_to_install) > 0:
    utils.install_packages(StrVector(names_to_install))

base = importr('base')
base.source('http://www.bioconductor.org/biocLite.R')
biocInstaller = importr('BiocInstaller')

biocInstaller.biocLite('SingleCellExperiment', suppressUpdates=True)
biocInstaller.biocLite('scran', suppressUpdates=True)
biocInstaller.biocLite('scater', suppressUpdates=True)

sceLib = importr("SingleCellExperiment")
scranLib = importr("scran")
scaterLib = importr("scater")
named_list = robjects.r.list(counts=counts)
sce = sceLib.SingleCellExperiment(named_list)

try:
    clusters = scranLib.quickCluster(sce)
    sce = scaterLib.computeSumFactors(sce, cluster=clusters, positive=True)
except rpy2.rinterface.RRuntimeError:
    print("fewer cells than the minimum cluster size, so no quickCluster is run")
    sce = scranLib.computeSumFactors(sce, positive=True)

----> sce = scaterLib.normalize(sce, logExprsOffset=1)

Everything runs fine except the last line which gives the error:

AttributeError: module 'scater' has no attribute 'normalize'

This is strange since that is what I am supposed to do actually according to docs.
I suspect that maybe I am loading a wrong scran version? I do not know though neither how to check the version number, nor how to install a specific one. Any suggestions would be greatly appreciated.
I tried the following:
sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)

And it is working. I am afraid though that it is a wrong way to go.
I tried to get the normalized matrix but it retrieved the initial raw counts one:
sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)
mtx_norm_r = robjects.r.assay(sce)
mtx_norm = np.array(mtx_norm_r)

So, something is definitely wrong here. Either I am retrieving the normalized matrix in the wrong way or the normalization was not run correctly. I tried the following method for retrieval but no luck:
mtx_norm_r = sceLib.normcounts(sce)

RRuntimeError: Error in assay(object, i = exprs_values) :
    'assay(, i="character", ...)' invalid subscript 'i'
  'normcounts' not in names(assays())

I found a way to get the version numbers using getNamespaceVersion() function:
robjects.r.getNamespaceVersion('scater')

The versions are:
scater - 1.8.4
scran - 1.8.4
SingleCellExperiment - 1.2.0

They are the newest as I understand.

I was able to get to the normalized values using logcounts function in the following way:
mtx_norm_r = sceLib.logcounts(sce)
mtx_norm = np.array(mtx_norm_r)

Since I used log_exprs_offset=1 the field that I needed to use is logcounts, otherwise it should be normcounts which I found here.

Maybe you know it but you can read the help page for each function using ?name_of_function on the terminal (There is no need to consult an online help page that might not match your version of the package).

But I can not check it in spyder

I don't know what is spyder but if it can run R is run you can.

It can not be done, or I do not know how anyway. It is rpy2 python library that translates objects between R and python and the documentation is often not clear and detailed enough, in my view, probably because not that many people use it. Spyder is just IDE for python from anaconda


text:
K mean clustering issue
So I have expression data set with 4 healthy sample and 4 disease, to determine the number of K I used mclust which i ran on the list of differentially expressed genes so I got around 7 cluster from Mclust which im using as an input for my K-mean .
Now my data set is like all the differentially expressed genes it contains both UP and down regulated genes so basically it may be UP in my healthy and Down in patient and vice-versa ,4 likely possibilities.
So now I have this doubt that from that I'm getting 7 cluster.Is that result is correct ?

What do you mean that you have 4 possibilites? Which constrasts did you do to calculate the differentially expressed genes? Did you compare healthy to disease? Then you only have 2 options.

There is generally no such thing as correct or incorrect when it comes to clustering. When you specify $K$, the program did it's statistical algorithms. Unless the the program had a bug, your results were statistically valid.
You will need to check the clustered labels visually. Only your biological knowledge can tell you the optimal $K$.

7 clusters may be correct in that it's the optimal number according to the Bayesian information criterion (or whatever value you had mclust optimize), but that doesn't make it necessarily the most sensible value. Have a look at your data and visually look at what might make sense. It's likely that the clusters are subdividing up and down regulated genes into various subgroups, which may or may not be useful, depending on your goals.

yes I used BIC ,and this is what my PI was asking if those cluster are even correct or make sense ...may be I would try some other clustering and see if get similar number I would let know...

As you're learning, clustering tends to be an iterative process.

very naive question when i define the number of k , the k mean clustering happens on my rows or columns or both ?

Rows are clustered using the values in the columns.


text:
Why is Clostridium difficile toxin A + B outdated in snomed?
I'm looking at the new version of the Snomed dictionary and the concept 707993000 Clostridium difficile toxin A + B (substance) is deactivated (and outdated) since version 20170731. Can anybody explain me, why this happened?
We have used this concept for coding of results of laboratory examinations in our software system. Is there an alternative Snomed concept usable for the purpose of coding of evidence of the Clostridium difficile toxin A + B ?


text:
using SNPs to identify mixed samples
Is there a way to identify mixed samples based on SNPs?
Example input (table of genotypes for multiple samples):
|      | sample1 | sample2 | sample3 |
|------|---------|---------|---------|
| rs1  | AA      | BB      | AB      |
| rs2  | AA      | BB      | AB      |
| rs3  | AA      | BB      | AB      |
| rs4  | AA      | BB      | AB      |

Desired output: Sample relationships/distances. In this example, sample1 and sample2 are distinct, but sample3 is likely a mix of two (biological or technical).
Based on this simple example, the calculation seems doable, but maybe it's impossible in more complex realistic scenarios.
There is a nice tool peddy that does something conceptually similar, but it expects pedigree information as well as genotype frequencies for both the input samples and population.

I updated the question. I didn't initially realize it was not clear.

Thanks, but this still remains a question that is mostly without context. What is the specific problem that you are trying to solve?
As presented here, those samples could be individuals and the problem is a phasing problem, rather than mixed groups. What platform is being used to carry out the genotyping? How far away from each other are the variants?

Yes. The general term for identifying subsets from a mixed group is deconvolution.
This is easier to do when there is haplotype-level information available (e.g. sequenced reads that span multiple SNPs), but it is possible to get some indication of sample mixes from unlinked SNPs alone.
Deconvolution is helped a lot by having proportional data in the input (e.g. frequencies of genotypes, rather than a single genotype). Most SNPchip platforms have multiple probes per variant, and can provide some indication of variant proportion when digging down into the raw data. Deconvolution is also a lot easier with some knowledge about the component mixes. It's much easier to simulate mixes than to guess at possible mixes.
Here's a paper about the application of deconvolution in forensics.
Here's a paper about a program called CIBERSORT which carries out a deconvolution of immune populations based on gene counts:
Here's another paper which discusses the idea in more general terms.
You might find more information on one of these Biostars posts:

Deconvolution Methods on RNA-Seq Data
Differential gene expession analysis in cell populations of mixed tumor and normal cells

I am familiar with the expression-based deconvolution, but that's a different problem. And you can have the same mix of cell populations regardless of how many individuals are mixed in a single library. Do you have any suggestions for SNP-based deconvolution?

Could you please clarify this in your question. It's not clear how much research you've done on this already, and whether or not you have haplotype-level information available.


text:
How to get the product of a CDS
I need the name of the protein in /product="protein_name" using bash commands. Beware, there is a lot of whitespace between lines.
FEATURES             Location/Qualifiers
     source          1..1266
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /isolate="Salem"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
                     /country="Australia: Taronga Zoo, Sydney"
     gene            1..1266
                     /gene="UROS"
                     /note="Derived by automated computational analysis using
                     gene prediction method: Gnomon. Supporting evidence
                     includes similarity to: 10 Proteins, and 100% coverage of
                     the annotated genomic feature by RNAseq alignments,
                     including 5 samples with support for all annotated
                     introns"
                     /db_xref="GeneID:100933602"
     CDS             60..908
                     /gene="UROS"
                     /codon_start=1
                     /product="uroporphyrinogen-III synthase isoform X1"
                     /protein_id="XP_003755055.2"
                     /db_xref="GeneID:100933602"
                     /translation="MRIPTSFSRGAWWKMKVLLLKDAKDDDSGLDPYLQELGSCGFEA
                     TLIPVLAFEFLSLPSFSEKLSHPEAYGGLVFTSPRAVEALQLSLERAGQLEVWKKSLK
                     EKWRLKPAYVVGQATASLVKKVGLAAEGEKCGNGEKLAEYICSRETPATLPLLFPCGT
                     LKGEVLPKTLKQKGILLESLTVYQKVQHPRMQDSLKSYFSTQGVPAGVTFFSPSGLTY
                     SLPHIQELSGASFDRIKFAAIGPTTARALASAGIPVSCTAEKPCPQDLAAGLQKALQP
                     SGPSSS"

Right now I have this:
while read p;
do cat "${p}.gb" | tr "\"" " "| egrep -v "/db_xref= taxon"| egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk '{print($2)}END{print("")}'
done <XM_codes.txt

Which results in this:
GeneID:100933602
UROS
uroporphyrinogen-111
XP_003755055.2

My output should be:
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

If you have GNU grep (default on Linux), you can do this with a single grep command using Perl Compatible Regular Expressions (PCREs):
$ grep -oP 'GeneID:\d+|/(gene|product|protein_id)="\K[^"]+' file | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

The -o flag tells grep to only print the matching portion of the line and the -P enables PCREs which let us use d+ for "one or more digits" and \K which means "discard whatever was matched before the \K".
Note that I'm not using cat here. The command from your question is what is known as a UUoC (useless use of cat) and is a surefire way of annoying any *nix geeks around ;).

If you don't have GNU grep, you can do something like this:
$ awk -F= '/GeneID:|(gene|product|protein_id)=/{print $2}' file |
    tr -d '"' | sort -u
GeneID:100933602
uroporphyrinogen-III synthase isoform X1
UROS
XP_003755055.2

Thank you so much! Also, what to do if the "uroporphyrinogen-III synthase isoform X1" string was on more than one line?

@AronvanBeelen ah, that is a very, very different situation. Is that even a valid genbank file? How would it look like? You would also need to remove the spaces on the second line, right? That is approaching the level of complexity where it is better to write a proper parser instead of hacking something together on the command line. You might want to ask a new question, showing an example that accurately represents your real data.

I've asked a new question, I hope you can help me! :)

If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t").  In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter.
 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

results:
GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2

The solutions here are more robust than mine since I rely on the structure of the file, but if you want a solution that is really easy to understand:
$grep -A 5 "CDS" file | grep -v "CDS" | grep -v "codon_start" | cut -f 2 -d \"
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2
GeneID:100933602

I take CDS line and 5 lines bellow, then I kick out CDS and codon_start (lines we are not interested in) and then I just bring what is between quote using cut. I also report the values in the order as they appear in the original file.


text:
How to correctly parallelise RSeQC scripts with GNU parallel?
I have a .bam, as ouput from STAR aligner,  from which I need to extract some info using RSeQC while using all the computational resources available to increase speed. To accomplish that I am using GNU parallel. Currently I am doing it like this:

Split the main .bam per chromosome, generating files with the suffix Aligned.sortedByCoord.out.REF
bamtools split -in Aligned.sortedByCoord.out.bam  -reference

Calculate the clipping profile
ls Aligned.sortedByCoord.out.REF_* | parallel --joblog clipping_profile.log --eta --bar --max-procs 20 --jobs 20 'clipping_profile.py -i {} -s "SE" -o RSeQC.clipping.profile/clipping_profile.{.}'

This will generate a .r, .pdfand a .xlxs for each chromosome, thus ending up with a lot of files. I am forced to split the .bam per chromosome as I am not able to find a --recstart nor a --recend in the compressed .bam.
Question
Is there a better way to split a .bam file using GNU parallel and RSeQC' programs (e.g. split_bam.py), such that one gets the outputs from a single input file?
Inspired by ole.tange's thread, I have also tried this, though unsuccessfully:
cat star_outputAligned.sortedByCoord.out.bam | parallel --block 1G -k --joblog splitbam3.log --max-procs 20 --jobs 20 --pipe 'split_bam.py -i {} -r ../../My_genomes/GCA.GRCh38/GRCh38_wgEncodeGencodeAttrsV28_all_rRNA.bed12 -o splitbam3'


text:
Find a cutoff value for genes that are expressed in single cell RNA-seq?
I want to find a cutoff value for each gene, above which we can consider a gene expressed.
The problem is that not all effectively non-expressed genes will have 0 counts due to sequencing errors for example.
My question is related to this one, but is concerned with single-cell data. The technical differences between bulk and single cell RNA-seq warrant a separate question for this problem. Also similar to this one, but focuses on gene expression and not differences between zero counts.

This cutoff would be applied to non-transformed values or to raw counts? What kind of counts do you have: at the transcript or at the gene level?

You have to be careful when setting these kinds of thresholds in single cell experiments. Using a mean or variance = cutoff may cause you to remove genes that define a rare subset of cells. For example, the genes derived from a subset of cells that represent 2% of your total population may have a low mean expression.

@Llopis I don't mind either way. The idea is to assign expressed / non expressed status to each gene (with some confidence). I have counts at the gene level.

@GWW Thanks for the comment. Yes I agree a simple approach won't do. Maybe for each gene we can look at the distribution of expression and choose the first minimum, or take a subset (cluster of similar cells, cell type) of cells and work with their mean / variance.

Not sure if this has anything to do with your question but I think the way how Monocle approach this question sort of makes sense to me. The normalization methods for Seurat and Monocle are somehow similar (correct if I am wrong, this is simply my personal experience), as @kohlkopf mentioned.
For monocle specifically, it has this function where you can display dispersion v.s. expression of the genes to decide the cut off. Here is an example of this plot that I generated when processing the PBMC data downloaded from the Seurat tutorial.
You can assign the cut-off manually but I doubt I answered your question on how to decide the cut-off value. Maybe you need to try a set the expression cut-off from a range to experiment what value may potentially answer a biologically relevant question in your context.
If you are interested in how monocle approach this, I think the Monocle introduction and Dave Tang's Blog are excellent in explain how it works.
Again, you are right that 0 could be sequencing error. But I would also suggest that, if you are looking at genes at very small expression level, scRNA-seq, at least droplet-based technique, could be really difficult in such context. I happened to hear someone trying to find lncRNA expression talking about how a "0" sometimes may be a true "0" due the low expression nature of lncRNA in scRNA-seq, etc.
Anyhow, not sure if I answer the question but I do hope it helps you somehow. I am not able to make comments but this is the best I can do.

Much of the Seurat workflow incorporates normalized, scaled, and mean-centered expression values.
Raw gene expression counts for a cell are normalized by total expression, multiplied by 10,000, and log transformed. Expression level for each gene is scaled by dividing the centered gene expression levels by their standard deviations, resulting in z-scores. It's not an ideal solution, but I find a threshold of 1 (in the realm of z-scores) to make sense. A z-score of 1 is one standard deviation from the mean, or what may be interpreted as baseline expression.
This is my 2¢. I would love to hear what other people are doing.


text:
Extract certain fields using from GenBank file using Bash script
I need to extract certain lines from a GenBank file on Linux.
Input file:
LOCUS       XM_012540068            1713 bp    mRNA    linear   MAM 11-JAN-2018
DEFINITION  PREDICTED: Sarcophilus harrisii cytochrome c oxidase assembly
            protein COX15 homolog (LOC100925550), mRNA.
ACCESSION   XM_012540068
VERSION     XM_012540068.2
DBLINK      BioProject: PRJNA167725
KEYWORDS    RefSeq; corrected model.
SOURCE      Sarcophilus harrisii (Tasmanian devil)
  ORGANISM  Sarcophilus harrisii
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Metatheria; Dasyuromorphia; Dasyuridae; Sarcophilus.
COMMENT     MODEL REFSEQ:  This record is predicted by automated computational
            analysis. This record is derived from a genomic sequence
            (NW_003816462.1) and transcript sequence (GEDN01126617.1) annotated
            using gene prediction method: Gnomon, supported by EST evidence.
            Also see:
                Documentation of NCBI's Annotation Process

            On Jan 11, 2018 this sequence version replaced XM_012540068.1.

            ##Genome-Annotation-Data-START##
            Annotation Provider         :: NCBI
            Annotation Status           :: Full annotation
            Annotation Name             :: Sarcophilus harrisii Annotation
                                           Release 102
            Annotation Version          :: 102
            Annotation Pipeline         :: NCBI eukaryotic genome annotation
                                           pipeline
            Annotation Software Version :: 8.0
            Annotation Method           :: Best-placed RefSeq; Gnomon
            Features Annotated          :: Gene; mRNA; CDS; ncRNA
            ##Genome-Annotation-Data-END##

            ##RefSeq-Attributes-START##
            assembly gap :: added 136 transcript bases to patch genome assembly
                            gap
            ##RefSeq-Attributes-END##
PRIMARY     REFSEQ_SPAN         PRIMARY_IDENTIFIER PRIMARY_SPAN        COMP
            1-136               GEDN01126617.1     1-136
            137-327             AEFK01003589.1     1071-1261
            328-450             AEFK01003589.1     2007-2129
            451-637             AEFK01003589.1     2869-3055
            638-805             AEFK01003589.1     6071-6238
            806-887             AEFK01003591.1     516-597
            888-1042            AEFK01003591.1     1431-1585
            1043-1156           AEFK01003591.1     4257-4370
            1157-1713           AEFK01003591.1     6191-6747
FEATURES             Location/Qualifiers
     source          1..1713
                     /organism="Sarcophilus harrisii"
                     /mol_type="mRNA"
                     /db_xref="taxon:9305"
                     /chromosome="1"
                     /map="unlocalized"
                     /sex="female"
     gene            1..1713
                     /gene="LOC100925550"
                     /note="The sequence of the model RefSeq transcript was
                     modified relative to its source genomic sequence to
                     represent the inferred CDS: added 136 bases not found in
                     genome assembly; Derived by automated computational
                     analysis using gene prediction method: Gnomon. Supporting
                     evidence includes similarity to: 10 ESTs, 11 Proteins, and
                     100% coverage of the annotated genomic feature by RNAseq
                     alignments, including 5 samples with support for all
                     annotated introns"
                     /db_xref="GeneID:100925550"
     CDS             53..1288
                     /gene="LOC100925550"
                     /note="The sequence of the model RefSeq protein was
                     modified relative to its source genomic sequence to
                     represent the inferred CDS: added 84 bases not found in
                     genome assembly"
                     /codon_start=1
                     /product="cytochrome c oxidase assembly protein COX15
                     homolog"
                     /protein_id="XP_012395522.2"
                     /db_xref="GeneID:100925550"
                     /translation="MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH
                     STLSEVALQSGRGAVSVPTKAAQRIMGRWLLVCSGTVAGAVILGGVTRLTESGLSMVD
                     WHLIKGMKPPTSQEEWEAEFQKYQQFPEFKILNHNMTLAEFKFIWYMEYSHRMWGRIV
                     GLAYVLPTVYFWRKGWLTHGMKGRVLALCGLVCFQGLLGWYMVKSGLEEKPDSYDIPR
                     VSQYRLAAHLGSALVLYCASLWTGLSLLLPRDKLPETRKLLQLRRFAHGTAGLVFLTA
                     LSGAFVAGLDAGLVYNSFPKMGDHWIPEDLFTFSPVLRNVFENPTMVQFDHRVLGIAT
                     VTAITGLYFFSRRIPLPRRAKMAVTTLLAVAYAQVALGISTLLTYVPTPLAAVHQSGS
                     LALLSIALWLMNELRRIPK"
ORIGIN
        1 gtggggcgac cgagggaggc cgctgcgagg ctgtgacttt gagactcgcg cgatgatgcg
       61 attgctgttt cctcctctga aagcggtgct gggcagtccg ggcgtccggc tgcaggctgc
      121 tgggggagct caccggtgcg gttgccatta tgggcaaaag ctctctctct tgcaaagaca
      181 gcacagtact ctgtctgaag tagctttgca gtctggaagg ggagcagtgt ctgtccccac
      241 aaaggctgcc cagcgcatta tgggtcgttg gctcttggtc tgcagcggaa cagtagcagg
      301 agcggtgatt cttggtggag tgactaggtt gacagagtct gggctctcca tggtagattg
      361 gcacttgata aaaggaatga aaccacccac aagccaagaa gaatgggaag cagagttcca
      421 aaaataccag cagtttccag aatttaaaat cttgaaccat aacatgacat tagcagaatt
      481 caagtttatc tggtacatgg aatactcaca cagaatgtgg ggtcggattg tgggcctagc
      541 atatgtcctt cccactgtct acttttggag aaagggctgg ctcacccatg gcatgaaagg
      601 tcgtgttctt gccctctgtg gtttagtctg tttccaggga ctgttggggt ggtatatggt
      661 gaagagtgga ttagaagaaa agccagattc ttatgacatc cctcgggtga gccagtatcg
      721 tctggccgct cacttgggct ctgcccttgt cctctactgt gccagcttgt ggactggact
      781 ctccctcttg cttcctcgag acaagttgcc agaaacccga aagctcttgc agctgaggcg
      841 atttgctcat gggacagcag gcctggtttt ccttacagct ctttcaggag cctttgtggc
      901 aggcctggat gctggacttg tatacaactc ctttcctaag atgggagacc actggatccc
      961 agaagatctc tttaccttct ctccagttct tagaaatgtc tttgagaatc caaccatggt
     1021 gcagtttgat catcgagtcc tgggaatcgc tacagtcaca gccattactg gactttactt
     1081 cttctccagg agaatccctc ttcctcgaag ggccaaaatg gcagtgacaa cattgctggc
     1141 tgtggcctat gcacaggtgg ccttgggcat tagtactcta ctaacgtatg ttccaactcc
     1201 tttggctgca gttcaccagt cagggtccct cgctttgctc agcattgccc tttggcttat
     1261 gaacgaactg aggcgaatcc caaagtaatg cccaggcagg tatttttctg gactcttccc
     1321 aggcttttga gagctctgta gagggcatga ggcttccctc ctgttgaatt aaagaattga
     1381 tgtggtagac agatcagaag atcaagaatg aagttactga tctgaggatt ctttttcaga
     1441 ctaaacaaac tactcttctt ccataagttg tcacttgctt tttaatgtgc tattaagatc
     1501 ccttccaact ttcactactc tgattctttt gatacagcaa taccactgat aggtctgtgt
     1561 ctcagagttg ggagggaggg aaggaaaaaa aaaacttttt gtacaaacgt acaaaactat
     1621 agtagttcct tttgtgggaa ctaaaataag gaaattttct ttctttttta tttattcgtt
     1681 ttcaacattt acttttataa gatttagatt ttc
//

My code:
tr "\"" " " < XM_012540068.2.gb | egrep -v "/db_xref= taxon" | awk -F= '/GeneID:|(gene|product|protein_id|translation)=/{print $2}' | sort -u

What I get:
 GeneID:100925550
 LOC100925550
 MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH
 XP_012395522.2
 cytochrome c oxidase assembly protein COX15

What I want:
 GeneID:100925550
 LOC100925550
 MMRLLFPPLKAVLGSPGVRLQAAGGAHRCGCHYGQKLSLLQRQH...(sequence goes on...)
 XP_012395522.2
 cytochrome c oxidase assembly protein COX15 homolog

Preferably, it has to be a mastercode for every GenBank-file. For a project I've got to loop through 78 GenBank-files, and extract the same data. I asked a similar question earlier, but the solution there doesn't work when the protein's name spans multiple lines.

Here is a script for achieving your result not so good but giving the desired output
#!/bin/bash
egrep -v "/db_xref= taxon"  $1| awk -F= '/GeneID:|(gene|protein_id)=/{print $2}' | sort -u > file1
sed -n '/product/,/protein_id/p' $1|grep -v protein_id |sed 's/.*=//' >file2
awk '/translation/,/ORIGIN/{if(/translation|ORIGIN/) next; print}' $1 >file3
head -3 file1 && cat file2 && printf "\"" && cat file3 |tr -d " \t\n\r" &&   echo && rm -rf file1 file2 file3

Save it in a file names as script1 and executing it as
`./script5 filename`

Gives following output
"GeneID:100925550"
"LOC100925550"
"XP_012395522.2"
"cytochrome c oxidase assembly protein COX15
                 homolog"
"STLSEVALQSGRGAVSVPTKAAQRIMGRWLLVCSGTVAGAVILGGVTRLTESGLSMVDWHLIKGMKPPTSQEEWEAEFQKYQQFPEFKILNHNMTLAEFKFIWYMEYSHRMWGRIVGLAYVLPTVYFWRKGWLTHGMKGRVLALCGLVCFQGLLGWYMVKSGLEEKPDSYDIPRVSQYRLAAHLGSALVLYCASLWTGLSLLLPRDKLPETRKLLQLRRFAHGTAGLVFLTALSGAFVAGLDAGLVYNSFPKMGDHWIPEDLFTFSPVLRNVFENPTMVQFDHRVLGIATVTAITGLYFFSRRIPLPRRAKMAVTTLLAVAYAQVALGISTLLTYVPTPLAAVHQSGSLALLSIALWLMNELRRIPK"

Hi Ammar, this script helps me a great deal! The only thing that I miss is the full string: "cytochrome c oxidase assembly protein COX15 homolog". Your script gives me: "cytochrome c oxidase assembly protein COX15"

@AronvanBeelen I have edited it again, hope it helps this time.


text:
Question: UCSC hub: autoscale option for multiwig track
I am making a hub on UCSC Genome Browser and I am aiming to display multiple bigWig files as a single track (multiwig container). I decided to use the autoscale option because I wanted a dynamic scale: y axis changing depending on the maximum values in the current observed window.
However, I am confused as to how to the relationship between:

autoscale on and viewLimits
autoscale on and type bigWig <#> <#>

For the "type" option, the official documentation mentions that: "The remote data files of type bigWig must declare the expected signal range for the data"
So my question is, are the following assumptions correct:

If I have autoscale on, I don't need to specify viewLimits
I still have to specify "type bigWig" but you no longer have to put a range (i.e. the <#> <#>)

If this is correct, my trackDb.txt file will look as such:
track data
container multiWig
configurable on
shortLabel label data
longLabel long label data
visibility full
priority 1
type bigWig
autoScale on
aggregate transparentOverlay
windowingFunction mean
showSubtrackColorOnUi on

    track sample1
    bigDataUrl path/to/sample1.bw
    shortLabel s_sample1
    longLabel s_sample1
    type bigWig
    visibility full
    alwaysZero on
    priority 120
    color 220,180,180
    parent data on

    track sample2
    bigDataUrl path/to/sample2.bw
    shortLabel s_sample2
    longLabel s_sample2
    type bigWig
    visibility full
    alwaysZero on
    priority 120
    color 220,180,180
    parent data on

    etc

Apologies, just updated my question, I hope it's clearer.


text:
How to run proteowizzard on linux using wine?
I am trying to run proteowizzard using wine to be able to convert raw mass spectrometry files to myML format. Following, these instructions` which seem terribly outdated. I am getting many issues. I am running Ubuntu 16.04.
We want to convert output from a ThermoFisher instrument to mzML or mzXML format.
> ~/software/proteowizzard>wine msconvert.exe
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:process:SetProcessShutdownParameters (00000380, 00000000):
> partial stub. err:module:import_dll Library VCOMP110.DLL (which is
> needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\baf2sql_c.dll") not
> found err:module:import_dll Library baf2sql_c.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library VCRUNTIME140.dll (which is needed
> by L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-convert-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-environment-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-filesystem-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library VCRUNTIME140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-filesystem-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library
> api-ms-win-crt-convert-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-environment-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library msparser.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:LdrInitializeThunk Main exe initialization for
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe" failed,
> status c0000135

Which issues did you find? You mention wine, so I assume you are using a linux computer to run a program for windows. Could you add more information about the requirements of this program?

It might be better to ask this on [ubuntu.se], but either way please [edit] and show us exactly what commands you ran and the specific errors you are getting. We can't help you if you don't tell us what the symptoms are.

I am wondering whether there is anyone here who runs proteowizzard on linux, who knows how to set it up.

Try this Docker container.
All the vendor formats work except T2D, YEP, and FID.

I sympathize ...
If it is a Mac just run 'PlayonMac' its a brilliant frontend for Wine without all the bugs and hassles it causes. Its free and available here


text:
What does "ALT+ cancer cells", the plus sign mean?
In many scientific papers, there is a plus sign:

... T-circles have also been detected in ALT+ (telomerase-negative) cells9 and ...

What does the "+" mean? Does it mean differentially expressed?

Here are a couple of papers (A and B) that explain ALT.

Telomere maintenance is a hallmark of cancer as it provides cancer cells with cellular immortality. A significant fraction of tumors uses the alternative lengthening of telomeres (ALT) pathway to elongate their telomeres and to gain an unlimited proliferation potential. Since the ALT pathway is unique to cancer cells, it represents a potentially valuable, currently unexploited target for anti-cancer therapies.
  B

The ALT+ denotes cell lines that are positive for the ALT pathway.


text:
Fastqc- Per Base Sequence Quality
I am trying to figure out how to interpret the "Per Base Sequence Quality"? What does Position in reads (bp) mean? Also in order to draw this box-plot graphics, more than one quality scores are needed so we can sort them and calculate the mean, median and quartiles. In this diagram, how these values were calculated for each position?

"in order to draw this box-plot graphics, ": what box plot graphics? "In this diagram": what diagram? Please [edit] your question and include an example of the data you are looking at. Ideally, also show us the command you used to generate it so we can reproduce it on our systems.

In this Page you can find good answer.

There is a pretty good description here. The position in the reads means literally the position in the reads, so position 2 means the quality of all bases on position 2 in the fastq file (all reads). In the description link, you can find how the boxplots and blue and red lines are calculated (and what they mean). They are all based on the total quality values of all bases on a certain position, all bases means all bases of all the reads on a particular position.


text:
How to format the Nexus file to create geographical maps with haplotypes in PopArt?
I'm trying to create a map in PopArt that allows me to visualise the distribution of several species. These are my Traits (=species) and GeoTags blocks:
BEGIN TRAITS;
    Dimensions NTRAITS=9;
    Format labels=no missing=? separator=Comma;
    TraitLabels SNE Lsp A B C D 3n 4n Lep;
Matrix
SNE1 1,0,0,0,0,0,0,0,0
SNE2 1,0,0,0,0,0,0,0,0
SNE3 1,0,0,0,0,0,0,0,0
SNE4 1,0,0,0,0,0,0,0,0
SNE5 1,0,0,0,0,0,0,0,0
SNE6 1,0,0,0,0,0,0,0,0
PUR_4n 0,0,0,0,0,0,0,1,0
GRA 0,1,0,0,0,0,0,0,0
OJE 0,1,0,0,0,0,0,0,0
ESC 0,0,0,0,0,0,0,0,1
LNM 0,0,0,0,0,0,0,0,1
PIN 0,0,0,0,0,0,0,0,1
URD 0,0,0,0,0,0,0,0,1
VIS 0,0,0,0,0,0,0,0,1
LYE 0,0,0,0,0,0,0,0,1
BEN_A 0,0,1,0,0,0,0,0,0
ALJ1_A 0,0,1,0,0,0,0,0,0
ALJ2_A 0,0,1,0,0,0,0,0,0
ESP2_B 0,0,0,1,0,0,0,0,0
ESP3_B 0,0,0,1,0,0,0,0,0
ALC_C 0,0,0,0,1,0,0,0,0
AVA1_3n 0,0,0,0,0,0,1,0,0
AVA2_3n 0,0,0,0,0,0,1,0,0
PER_3n 0,0,0,0,0,0,1,0,0
TRI_3n 0,0,0,0,0,0,1,0,0
VDI_3n 0,0,0,0,0,0,1,0,0
VIL_3n 0,0,0,0,0,0,1,0,0
CVE1_4n 0,0,0,0,0,0,0,1,0
CVE2_4n 0,0,0,0,0,0,0,1,0
LNA1_4n 0,0,0,0,0,0,0,1,0
LNA2_4n 0,0,0,0,0,0,0,1,0
PRA1_4n 0,0,0,0,0,0,0,1,0
PRA2_4n 0,0,0,0,0,0,0,1,0
VEN1_4n 0,0,0,0,0,0,0,1,0
VEN2_4n 0,0,0,0,0,0,0,1,0
ZAR_4n 0,0,0,0,0,0,0,1,0
SFI1 0,0,0,0,0,0,0,0,1
SFI2 0,0,0,0,0,0,0,0,1
SIB 0,0,0,0,0,0,0,0,1
TIS1_D 0,0,0,0,0,1,0,0,0
TIS2_D 0,0,0,0,0,1,0,0,0
TIS3_D 0,0,0,0,0,1,0,0,0
TIS4_D 0,0,0,0,0,1,0,0,0
TVI_D 0,0,0,0,0,1,0,0,0
PVRFAS0001 1,0,0,0,0,0,0,0,0
PVRFAS0002 1,0,0,0,0,0,0,0,0
PVRFAS0024 1,0,0,0,0,0,0,0,0
PVRFAS0025 1,0,0,0,0,0,0,0,0
PVRFAS0027 1,0,0,0,0,0,0,0,0
PVRFAS0028 1,0,0,0,0,0,0,0,0
;
END;
Begin GeoTags;
    Dimensions NClusts=29;
    Format labels=no separator=Spaces;
    [ClustLatitude 37.15 37.77 37.13 39.37 37.11 40.60 40.59 37.04 40.61 38.71 36.70 43.48 40.18 40.24 40.57 40.26 37.10 37.86 38.52 36.56 40.23 38.01 39.58 38.75 40.54 40.55 38.70 40.04 40.69;
    ClustLongitude -3.70 -3.02 -3.52 -3.72 -3.03 -4.33 -4.33 -3.43 -2.44 0.05 -5.35 3.55 -6.55 -2.06 -4.16 0.26 -3.26 -1.56 -3.56 -4.86 -2.07 -2.87 -3.87 0.18 -4.18 -4.19 -0.49 -7.49 -2.59;]
Matrix
SNE1 37.044888 -3.135564
SNE2 37.044888 -3.135564
SNE3 37.044888 -3.135564
SNE4 37.044888 -3.135564
SNE5 37.044888 -3.135564
SNE6 37.044888 -3.135564
PUR_4n 37.133158 -3.524806
GRA 36.708272 -5.348794
OJE 36.562847 -4.853522
ESC 40.5659132 -4.2236186
LNM 40.609845 -4.344144
PIN 38.752157 0.180615
URD 39.370381 -3.720659
VIS 38.520059 -3.564201
LYE 39.580131 -3.870267
BEN_A 38.712603 0.052372
ALJ1_A 38.704010 -0.490801
ALJ2_A 38.704010 -0.490801
ESP2_B 37.865615 -1.569499
ESP3_B 37.865615 -1.569499
ALC_C 40.264366 0.267744
AVA1_3n 43.490194 3.563363
AVA2_3n 43.490194 3.563363
PER_3n 40.610461 -2.447806
TRI_3n 40.697087 -2.595488
VDI_3n 40.244599 -2.067766
VIL_3n 40.231843 -2.078962
CVE1_4n 40.559003 -4.199853
CVE2_4n 40.559003 -4.199853
LNA1_4n 40.597289 -4.330246
LNA2_4n 40.597289 -4.330246
PRA1_4n 37.113314 -3.035071
PRA2_4n 37.113314 -3.035071
VEN1_4n 37.155124 -3.703330
VEN2_4n 37.155124 -3.703330
ZAR_4n 40.547662 -4.184395
SFI1 40.042132 -7.498515
SFI2 40.042132 -7.498515
SIB 40.182598 -6.551259
TIS1_D 37.771593 -3.023286
TIS2_D 37.771593 -3.023286
TIS3_D 37.771593 -3.023286
TIS4_D 37.771593 -3.023286
TVI_D 38.0140003 -2.8751617
PVRFAS0001 37.100000 -3.260000
PVRFAS0002 37.100000 -3.260000
PVRFAS0024 37.100000 -3.260000
PVRFAS0025 37.100000 -3.260000
PVRFAS0027 37.100000 -3.260000
PVRFAS0028 37.100000 -3.260000
;
End;

If I upload my Nexus file like this, the program throws an error saying "Too many rows in GeoTags block." However, I have 50 samples in both Traits and GeoTags. The same happens if I uncomment (in Nexus format, comments are in brackets) the ClustLatitude and ClustLongitude files.
If I set the "labels" flag on "yes", it takes the traits as the geographical locations and tries to find locations in the world starting with the trait names.
If I comment the "GeoTags" block, PopArt correctly outputs a network.

What am I missing? Does somebody have a Nexus example where the traits do not correspond to geographical location of the samples?

I don't know this format but you seem to have two BEGIN and two ENDs. To me it looks like you have two datasets in one. Also one is in caps and the other is not. Are you sure that the file is in the right format? Did you check it (if there is any tool to do that) ?

This program uses a custom version of the Nexus format. Here is the manual: http://popart.otago.ac.nz/doc/popart.pdf and it uses exactly the "begin" and "end" blocks I'm using - weird format, I know, but I'm trying to stick as much as possible to the example given...

@Llopis Additional info here (Supplementary data of the PopArt paper) and here (Nexus example in the PopArt site)

In the pdf you link in the first comment the traits data is formatted different (see the corresponding section of "Traits Tables"). Are you sure you should provide all these data in a single file? The error seems to imply that when an END is found it should be the end of the file or something alike.

Ah you mean the indent? That's a copy-paste error, I'll edit that. In the original file I indented the block between "BEGIN TRAITS" and "Matrix". Note that if I comment the GeoTags block, a network is produced correctly. The manual notes: " I've separated the individual blocks with text describing the important points, but they should all be included in a single text file, as shown in here."

No, I mean that the section "traits Tables" describe a different format (coma separated values I think) than the one you have, but trait information is described in two formats.

authors: Обедающий философ


text:
Random addition method for phylogenetic tree reconstruction
I have been working on nucleotide and protein data to reconstruct maximum parsimony trees using MEGA software to get parsimony tree, via "random addition" method is used to generate the initial tree(s).
I need to know that when randomly chosen sequences are added to reconstruct an initial tree, on what basis the best tree is chosen??
I would appreciate if anyone can guide me about this method or suggest me some good literature about this in order to make me understand the concept of this method.

Is there something wrong with question?

Nothing is wrong with this question, but people do this on their own time and are not working to solve other's problems, and some might not be on the same time-zone as you. Did you try reading the original article ?

Yes, I have gone through it but it did'nt help me solving my question.

MEGA manual may help to explain the criteria used in selecting the final maximum parsimony tree.
The random tree additions refers to the generation of random initial trees for a heuristic search to find MP trees. In this case, a tree is generated by randomly selecting a sequence and adding it to the growing tree on a randomly-selected branch.
Instead of using the random tree additions you can use the min-mini algorithm for the initial tree set.

Thankyou for the response, but the thing is that I need to understand how can an initial tree be constructed to be used for heuristic searches to find parsimonious tree?


text:
Boxplot of kmean cluster
 gene   H4       H3       H2       H1       L8       L7       L6       L5       L4       L3       L2       L1          clust
AAGAB  7.540712 7.296325 7.030218 7.408974 8.269541 7.758752 9.253013 8.611002 8.705045 8.940146 8.855133 8.488776     3
ABCA1  7.112951 6.213973 7.200973 6.784828 8.044299 9.832347 8.558342 7.134570 8.993938 5.794208 8.271960 9.179501     5
ABCA6  1.822567 2.209167 3.379618 2.813433 3.180290 3.741516 3.060194 2.078485 6.536346 6.881403 5.927874 5.278854     2

This is my output after doing kmean clustering so at the end I have my cluster number assigned to each gene. Now I want to show how the expression of my sample varies across each cluster, I am not sure how to transform the data and do it in ggplot2 as of now i doing something like this in base R plot
boxplot(da$H4 ~ k$cluster,
        xlab='Cluster', ylab='Genes',
        main=' Genes by Cluster')

For each sample I am checking, is there any way to show that with all the samples together, any suggestion or help would be highly appreciated.

Assuming da is a data-frame with the cluster as factor you can do the following:
ggplot(da) +
  geom_boxplot(aes(cluster, H4)) +
  xlab("Cluster") +
  ylab("Genes") +
  ggtitle("Genes by Cluster")

okay i did something similar  made factor of the clust column then melt  using  melt(da, id.vars=c("clust")) the plot , I got all the cluster expression together..


text:
How to query the Human Microbiome Project (HMP) to find all subjects with both 16s and WGS workups?
I am looking for a query to run on the HMP database that will return all subjects who have had BOTH 16s and whole genomes sequence (WGS) workups.  I am currently using this query...
file.matrix_type = 16s_community OR file.matrix_type = wgs_community AND (file.format = FASTA OR file.format = FASTQ)

and sifting through the results by adding them to my shopping cart, browsing the file types until I have added all of the samples from one subject, and then I check to see if I have at least one 16s and one WGS.  So far, I have not found one person who has had both 16s and WGS where both files were publicly accessible, and most do not have WGS at all.  There must be a way to get subjects with both workups done to make this faster.

After trying a few different queries, I don't believe there is a more specific one than this and I found that by searching and clicking through all of the results given by this query that eventually I found patients that fit the criteria of having both WGS and 16s.