text:
Making PLINK compatible files from VCF file without phenotype information
I have a big VCF file that I need to convert to, preferably, bed/bim/fam files that are readable by plink. Currently using plink 2. I am aware that this version of plink can be used to convert VCFs into binary peds using
plink --vcf file.vcf.gz --make-bed --out out
Or something similar. I also know that in this case, plink will automatically fill the phenotype and sex columns in the bed and fam files with zeros. I know there is a --pheno flag that allows specifying a txt file for the phenotypes and a --make-pheno flag that can be used for the same effect.
How exactly does this work? How can I encode the phenotype and sex information correctly? Can I just use --make-bed and then manually make a phenotype/sex file to use down the line with the --pheno flag? Or would I do something like
plink --vcf file --pheno phenotype_file --make_bed --out out_file
directly? In this case, will there not be a problem with the fact that the generated .bed in the first step does not have this information? What about with sex?
Another option I can think of is just edit the generated .fam to include the correct information and use it downstream, but that information will still not be in the .bed file.
I have no family relationships, so the FID fields would just be 0s (I guess I can add the --const-fid flag for that effect when I'm converting, right?).

cross-posted https://www.biostars.org/p/491146/

You can use --pheno to specify a phenotype file you want to use, the 1st and 2nd columns being FID and IID. These are just used on the fly and cannot be written to the BED files as I understand. --make-pheno is used to write binary phenotypes, which is unsuitable if your phenotype is actually continuous.
But I still do not quite understand what you want to do here.

I was thinking about the .ped file, not the .bed, hence the confusion. The .bed contains only information about the SNPs (unlike the .ped). Everything working fine now.

Sure, by BED files I mean the PLINK BED format, which includes the three files types, not the .bed file.


text:
extract chimeric and multimap reads from bam file
I am trying to extract all chimeric and multi-map reads from either SAM/BAM file. Is there any simple command to do that? Can I use htslib for parsing sam/bam files and extracting based on any flags. Any help on how to start will be appreciated.

This is how I do it.
samtools view -@5 -f 0x800 -hb /path/sample.bam > sup.bam

Then if you want it as a fasta
samtools view -F 0x1 -hb sup.bam | samtools fasta -F 0x1 - > sup.fa

Thanks !! It works great. I will use samtools source code to write a small program to extract the reads based on flag.

Pysam wraps it pretty well, just fyi

source: Stack Exchange

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

sort: votes


text:
how to get access to genotypes and phenotypes used for a GWAS

I'm a master's student working on genomic prediction of complex traits using deep learning.
i'm looking for a dataset of human genotypes and phenotypes that has been used for a GWAS.
The only thing i can get from the gwas i found is the summary statistics not the raw data.
Does anyone know if it's possible to get such data (for free ) for a research project?

Thanks.

How about:

dbGaP - https://www.ncbi.nlm.nih.gov/gap/ - in case you are searching for supporting datasets I leave its entry in BiŌkeanós

https://www.ddbj.nig.ac.jp/jga/ - the access might be a bit more complicated.

https://www.internationalgenome.org, not particularly large, in their FAQ they state that the individuals sequenced have not self-reported any diseases - perhaps you might want to add those with to your negative set.  Again, its entry in BiŌkeanós.

If you are looking for more comprehensive information, you could try biobanks - e.g. https://www.ukbiobank.ac.uk. These are usually open (or reasonably priced) for research but the approval process can take time.
A side note - some of GWAS catalogues will have a link to the study where the association originated from, e.g. here. You can sometimes find the datasource there.

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


text:
Reliably create female and male individuals in a -ped file for plink
For a study I need some kind of "decoy" GWAS-data. For this I created a panel with some genes, the same ones as in the real data. From NCBI I fetch some of their SNPs and create my own artificial .ped and .map file. Furthermore I save for each snp if it lies on the X-Chromosome.
I then go on and randomly create female individuals by giving them heterozygous alleles on their X-Chromosome genes. Males only get homozygous ones.
But when I then run pLink with the --check-sex argument I get a sexcheck file where there's only zeroes in the SNPSEX column. So do you know how --check-sex works that I can then create decoy-individuals that reliably (as much as possible) get detected as a given gender?

I was able to create a sample map/ped files that result in SNPSEX column being populated in plink 1.19 using 1000genomes Xchrom - found in their FAQ: https://www.internationalgenome.org/faq/where-are-snps-xymitochondrial-chr/.
The relevant conversion command ./plink  --check-sex --vcf ../ALL.chrX.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.vcf --maf 0.05 --out ped_map_chrX/chrmX --recode .
Even in this dataset, SNPSEX is 0 (ambiguous) in 400 (~16%) of cases. One could start pruning from here to arrive at a suitable minimal set of SNPs that are still recognised in plink.

Ok, thanks! I'll try that out.


text:
How to Interpret Author's Methods in this article?
I have to reproduce a result for school reasons (bioinformatics exam) and I was assigned a certain article. In particular I have to reproduce what some authors did in: doi: 10.1371 / journal.pone.0217889. From what I understand, they make functional enrichments in respect of three variants of protein disorder. They consider vectors whose components are proteins (belonging to a certain type of disorder flavour) in a functional class divided by the total number of proteins in that functional variant. The problem is that if I go for an enrichment in GO terms for a certain set of proteins, the same GO terms do not necessarily come out for another set of proteins. So my question is: in relation to the paragraph 'Clustering of the variants of disorder through functional profiles' how do I construct the vectors that they then use in clustering?
Thanks in advance
The full quote to the article is
Deiana A, Forcelloni S, Porrello A, Giansanti A. Intrinsically disordered proteins and structured proteins with intrinsically disordered regions have different functional roles in the cell. PLoS One. 2019; 14 (8): e0217889. Published 2019 Aug 19. doi: 10.1371 / journal.pone.0217889

I feel like it's for a given GO term, count the number of proteins belonging to a protein variant group in that GO term and divide it by the total number of proteins belonging to that variant.

Yep:up to there no problem. The problem is that when one feeds sets of different proteins, the GO terms are not comparable. I had an idea: at this point it would be better to take the entire proteome and see how it is noted and then divide by protein variants. At that point, GO labels would be obtained with which to build three vectors of GO terms. I'm doing so then I don't know how close it comes to the strategy adopted by the authors also because PANTHER is a bit dated DB and to be honest I don't like it very much. I think I will use g: profiler


text:
Explanation for RNA-seq samples not clustering in PCA as expected
A colleague is analysing RNA-seq data - the study design is 2 treatments, 3 replicates, 3 tissues. In their PCA plot the samples clustered neatly by tissue. Except for two samples - two tissue samples originating from the same animal clustered with the wrong tissue.
I suggested that the two samples had likely been switched accidently somewhere along the way. They think its explained simply by biological variation which I'm not convinced of. Are there any technical explanations for why samples could cluster on a PCA in such a way?

If the two odd samples each belong to the other's cluster, swapped labels seems like a reasonable hypothesis. You could look for tissue-specific genes (using something like GTex) to see if the expression patterns match the purported tissue as well. If you find high expression of liver-specific genes in a sample labelled as brain tissue and high expression of brain-specific genes in the liver sample, for example, you have good evidence to support the label swap hypothesis.

I agree with @NuclearHoagie, you can look at https://bgee.org/, they have curated per-tissue expression data of many species (if you have something at least a bit related to any of the models)

Biological variation is not going to make one tissue look like another, unless the two tissues are incredibly close, like adjacent parts of the brain.  I bet your PC1 determines a really hefty proportion of the variability, like 80%?  Biological variability isn't going to be nearly that strong.

It's hard to say what is going on without looking at the plots. But it is a reasonable guess that samples are mislabelled and I think it is worth checking the experimental log carefully. It is also possible that the two samples are of low quality. People often use PCA to check and remove outliers from their RNA-seq libraries. So you may also check the RNA-seq library quality for those two samples.

Check for mislabelled samples

For the mismatch you can check if there is a marker gene for these 3 samples for that particular cohort/tissue that you actually think should be expressed in this tissue only and not in other cohort.

Check for outliers

To check the logs of mapping reads of the samples that aren't clustering within expected tissue it might be due to low reads. You can use hierarchical clustering and project it in the form of a heatmap. Check if these two samples have varying expression levels to its corresponding biological replicates.

Have you corrected for batch effects? This is usually a reason, as it can mean that the variation being observed from your PCA is due to other factors rather than tissue and/or treatment as you expect - of course, another can be they are incorrectly labeled !
If the samples were sequenced at different times, or by different services, this can often happen. If that is the case, what you will usually find are samples being clustered based on the sequencing service, for example.
Use ComBat or limma packages in R to fix variations in the data that may be due to this, and rerun the PCA. Best of luck.


text:
Filtering Sequences (entries) by headers ID from a Fasta file database
First of all, sorry if this question has been posted previously (I could not found a solution accessing the previous Q&A).
I have a fasta file as follow:
>LNIV02000036.519060.520603 Pseudomonas aeruginosa
GAACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAA...

>LDLH01000045.361.1876 Pseudomonas aeruginosa
AGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATG...

>LNMQ01000087.5825.7394 Enterococcus faecium
TTTTTATGAGAGTTTGATCCTGGCTCAGGACGAACGCTGGCGGCGTGC...

...

I would like to filter sequences or extract entries for a specific species, let's say Pseudomonas (I have others in mind...). I tried a few things but I think the main problem is my skills associated with bioinformatics. I would like to create a file content Pseudomonas entries (~9k entries in the current database). Any suggestion or tip?
Thank you for your time.
Best,
Fabricio

You can do this using seqkit as follows:
seqkit grep -r -n -p '.*Pseudomonas.*' temp.fa

To explain a little, seqkit grep will allow you to search FASTA/Q files by sequence name or sequence itself. In this instance:

-r tells that the pattern is a regular expression
-n to match by full name instead of just id
-p to specify the regular expression pattern to search

There are a bunch of other tools that come with seqkit that can prove indispensable for basic sequence analyses.

Hi @vkkodali, thank you for your time posting it. I did not know about seqkit (living and learning). I'm testing it out and looks like a great tool. Thanks again.

Hi again @vkkodali, another question (sorry, I'm a novice in the world of bioinformatics). Is there a way to retrieve the whole sequence header or ID using seqkit? I filtered the sequences that belong to Pseudomonas and the fasta file contains 38K entries of Pseudomonas. I would like to extract the sequence ID to count how many species of Pseudomonas are in this file (38K entries).
illustration:

AB001446.1.1538 Pseudomonas amygdali pv. mori
AACTGAAGAGTT...

AB001440.1.1538 Pseudomonas coronafaciens pv. atropurpurea
AACTGAAGAGTTT...
...

Thank you for your time and attention.

So you just want the AB001446.1.1538 portion of the header? So many ways to do that... Easiest would be something like grep '^>' in.fa | cut -f1 -d ' ' | sed 's/>//'. Another way is to do seqkit faidx in.fa and that creates a in.fa.fai which is essentially a tab-delimited file with the first column containing the identifiers.

Hi, can be the whole header or portion of it. The portion is the genus and species name (e.g., Pseudomonas putid). I would like to extract the species name from the filtered Psedumonas file (e.g., Pseudomonas amygdali pv. mori, Pseudomonas coronafaciens pv. atropurpurea etc.). I would like to know how many different species of Pseudomonas are there. Thanks again for your time. Fabricio

This can be done using a combination of grep and sed but I'd be making assumptions about how the organism names in the headers for all your 9k sequences are formatted. Perhaps you should describe the entire issue with an example of what you are starting with and what you would like to see in the end and I can update my answer?

Hi @vkkodali, due to the limitation of characters in the comment section here. I wrote an answer which actually is a follow-up (explanation) for the same issue. Thank you for your time. Fabricio.

Using Biopython:
from Bio import SeqIO

for record in SeqIO.parse("input.fa", "fasta"):
    if record.description.split()[-2:] == ["Pseudomonas", "aeruginosa"]:
        print(record.format("fasta"))

Hi Chris, right after the post I realize that the issue was not properly described. So, this database has ~9k entries of Pseudomonas species, i.e., hundreds of P. species. The command above asks for the name of each Pseudomonas entry, correct ("Pseudomonas", "X")? Is there a way to retrieve the Pseudomonas entries without specifying all Pseudomonas names (I do not have a list name)? Thanks.

@Fabricio Hi Fabricio, do you understand the above code? If so you can easily modify the if line to search for entries with "Pseudomonas"

@Fabricio Add this at the top of the code, after the other import statement: import re . Change the if expression to this: if re.findall(r'pseudomonas', record.description, re.I): , in order to search in the fasta header for the word pseudomonas, in case-insensitive manner.

"describing the entire issue".
I have a database (silva 16S formatted by DADA2 - fasta format) that contains ~300k entries for bacteria (16S rRNA gene). e.g.,

AB001441.1.1538 Pseudomonas syringae pv. broussonetiae
AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGG...
KR054096.1.1224 Marinococcus halophilus
GAAGGTAAAGGCTTCCCCAAGGGCGACGATGGGTAGCCGACTTGA...
AB681405.1.1478 Paenibacillus barengoltzii
GACGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAGCGGAGT...
...

I would like to retrieve sequences for some genus, e.g., Pseudomonas. After running the seqkit command (mentioned by you above) I ended up with ~9 of Pseudomonas entries in a fasta file (example below).

AB001441.1.1538 Pseudomonas syringae pv. broussonetiae
AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA...
AB001442.1.1538 Pseudomonas amygdali pv. eriobotryae
AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA...
...

This pool of 9k sequences of Pseudomonas has several species of Pseudomonas (e.g., P. putida, P. syringe Et Cetera). Since some species of Pseudomonas have multiple sequences in the database (e.g., P. putida has 1.8k entries) I would like to know how many different species of Pseudomonas are there (in the 9k Pseudomonas entries), by names (i.e., in the filtered Psedumonas fasta file). At the end of the day, I would like to say that in this database I have 9k entries of Pseudomonas which represents 'X' number of Pseudomonas species.
What I'm thinking... Retrieve the full ID (AB001442.1.1538 Pseudomonas amygdali pv. eriobotryaeall) or if possible the portion of the ID (headers) that contains the species names (e.g. Psudomonas amygdali pv. eriobotryaeall) and with this, I can list the species names and count them. Does this explanation help to understand what I'm trying to do?
Let me know if you need more info.
Thank you for your time.
Fabricio


text:
Proper scatter/gather of bam files
I'm trying to do some calculus in parallel with BAM files. For that, I need to split the input bam file in chunks, process them in parallel and join them at the end. But I'm facing a read duplication issue on the split points.
I'm trying to split the input files by equally sized intervals (computed by GATK SplitIntervals), but that means the one of the files will contain all reads that covers the last base of the interval and the next file will contain all the reads that covers the following position. This means that we are putting in both files all the reads that crosses the boundary.
When we have to merge the results back into one file (via merge, sort or GatherBamFiles), all this repeated reads go into the result file, givin us a bam file with some extra reads that happen to have the same ID as others (because, in fact, they are the same).
Which would be the way to go to accomplish this task? How could I gather the files without duplicates? Or how could I split the files without putting the same read in more than one file?
For reference, these are the commands I'm trying to check that I'm duplicating reads. I extract a one base interval from a file, the next base interval from the same file and finally counting the number of reads in the files:
$ samtools view -c test.bam 1:18454-18454  # Reads covering position 18454
251
$ samtools view -c test.bam 1:18455-18455  # Reads covering position 18455
250
$ samtools view -c test.bam 1:18454-18455  # Reads covering 18454-18455 range
251

$ samtools view -b -o test.1_18454_18454.bam test.bam 1:18454-18454
$ samtools view -b -o test.1_18455_18455.bam test.bam 1:18455-18455

$ samtools view -c test.1_18454_18454.bam  # Reads extracted that covers position 18454
251
$ samtools view -c test.1_18455_18455.bam  # Reads extracted that covers position 18455
250

$ samtools merge test.1_18454_18455.merged.bam test.1_18454_18454.bam test.1_18455_18455.bam
$ samtools view -c test.1_18454_18455.merged.bam  # Reads merged that covers 18454-18455 range
501


text:
Nanopore data clarification regarding merging of samples
We had same set of sample ran on nanopore platform twice.
The issue on hand is.
First run was done it was nearly done but the power went off.
Second run was done completely with enough sequencing depth.
Now the question I ran the analysis with the second set which was complete that set of result I have. But so far I have't ran the first part.
Question1 : if I run the first part do i get major differences? in the result
Question2 : would it be advisable to merge both run since they are now can be called as twi different batches plus would be bring drastic changes in the output ? And in terms of tool for merging I googled which shows poretools.

It is entirely unclear which type of analysis you aim to do. You probably should be more specific when asking questions.
That said, given that the run was going perfectly fine until a power issue stopped it, all data produced should be fine and you can indeed combine those runs.
You probably shouldn't bother with poretools. If you are using fastq files to start your analysis you can just cat the basecalled files together. Also after alignment, you can merge the bam files. If your analysis requires the fast5 files then you can probably just put these in the same folder, but you should tell us which analysis you are doing.

So we have gut mircobiome sample ..between two groups. For classification Im using kraken2 and the output is used for downstream analysis using phyloseq


text:
Global alignment between two sequence X and Y with maximum number of identical matches
If 2 protein sequence X of length m and Y of length n and if there are several highest scoring global alignments, then I want to get one alignment that has largest number of columns in which a letter from X is aligned to an identical letter from Y. I need an algorithm to find alignment with maximum number of identical matches. Can anyone help me?

There are a lot of possible ways to count the number of perfect matches in an alignment, and the details will depend on which alignment tool you are using.
But, if you use the built-in pairwise aligner in BioPython (Bio.Align.PairwiseAligner(), https://biopython.org/docs/1.75/api/Bio.Align.html), then the easiest way is to use the formatted string representation of the alignment, which uses the pipe/vertical-bar character ('|') to indicate perfect matches.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
alignments = aligner.align("GAACT", "GAT")
max_matches = 0
for i, a in enumerate(alignments):
    align_str = f'{a}'
    align_str = align_str.split('\n')
    match_str = align_str[1]
    match_count = match_str.count('|')
    mismatch_count = match_str.count('.')
    print(f'Number of matches: {match_count}')
    print(f'Number of mismatches: {mismatch_count}')
    if match_count>max_matches:
        max_matches = match_count
        max_ind = i

best_alignment = alignments[max_ind]

In this example, both alignments have the same number of perfect patches (3), but that may not always be the case, depending on the input sequences and the scoring parameters used. The example just uses the default scoring parameters, but in general, you probably want to choose parameters appropriate for your type of alignment.


text:
STRINGdb: reverse map STRING identifiers to gene names?
Is there a functionality within STRINGdb package or in another environment to reverse map the STRING identifiers to Gene names?
For example, the entire STRING network can be downloaded using the commands:
library(STRINGdb)
string_db <- STRINGdb$new(version="11", species=9606,
                          score_threshold=0, input_directory="" )
full.graph <- string_db$load_all()
netFull
# IGRAPH c07d6e9 UN-- 19354 11759454 --
# + attr: name (v/c), combined_score (e/n)
# + edges from c07d6e9 (vertex names):
#  [1] 9606.ENSP00000000233--9606.ENSP00000272298 9606.ENSP00000000233--9606.ENSP00000253401
#  [3] 9606.ENSP00000000233--9606.ENSP00000401445 9606.ENSP00000000233--9606.ENSP00000418915
#  [5] 9606.ENSP00000000233--9606.ENSP00000327801 9606.ENSP00000000233--9606.ENSP00000466298
#  [7] 9606.ENSP00000000233--9606.ENSP00000232564 9606.ENSP00000000233--9606.ENSP00000393379
#  [9] 9606.ENSP00000000233--9606.ENSP00000371253 9606.ENSP00000000233--9606.ENSP00000373713
# [11] 9606.ENSP00000000233--9606.ENSP00000248929 9606.ENSP00000000233--9606.ENSP00000356737
# [13] 9606.ENSP00000000233--9606.ENSP00000338766 9606.ENSP00000000233--9606.ENSP00000386135
# [15] 9606.ENSP00000000233--9606.ENSP00000319192 9606.ENSP00000000233--9606.ENSP00000432568
# + ... omitted several edges

V(netFull)$name will provide all the STRING identifiers within this network netFull. How to get the corresponding gene names?
For example, string_db$map() maps the gene identifiers of the input dataframe to STRING identifiers. It returns the input dataframe with the "STRING_id" as an additional column. Now, I would like to explore in the reverse direction going from STRING identifiers to gene names.

I don't know if it's too late, but I met same question and got an answer.
You can use add_proteins_description(). It works in same way with map(). If you put a data frame having STRING_id column, it adds preferred_name, protein_size and annotation. You may find your answer on preferred_name column.
require(STRINGdb)
string_db <- STRINGdb$new()
string_db$add_proteins_description(data.frame(STRING_id = [vector containing string ids]))

This answer does not rely on STRINGdb package.
STRING provides supplementary identifier mapping files here: https://string-db.org/mapping_files/.
It looks like you are interested in human STRING name one - https://string-db.org/mapping_files/STRING_display_names/human.name_2_string.tsv.gz. You could read it as a separate lookup.


text:
About optimal distance (or dissimilarity) measure in proportion data
My data is proportional data, but variables are correlated with each other. So I am wondering what kinds of distance or dissimilarity we should use to take correlation into account.
If we would like to use data with proportion (e.g.,total sum scaling or compositional data), what kinds of distance or dissimilarity can we use?
Bray-Curtis is one of them, but does this make sense to use arcisin or logit transformation and use Pearson correlation?
In continuous variables, we can use Pearson and Spearman distance. In our case, what kinds of distance can we use? Pearson distance cannot be applied to data with 0-1 scale. Logit transformation -> Pearson?

I would suggest reading a little on the topic; my understanding is that distance measures can be pretty robust to proportion/compositional data. Here is a relevant textbook chapter:

Other distance measures assume that the data are proportions ranging from zero to one, inclusive.

Among these are indeed Bray-Curtis as well as Jaccard distance (and others). I think that you will want to use one of the more traditional distances rather than a Pearson/Spearman, though they are getting at the same problem. A transformation might help, I suppose, but I would prefer to choose a simple but predictable distance measure rather than a correction on a possibly unsuitable distance that might have unpredictable behavior.
It's a little hard to tell from your question, but it seems that you are wanting to in some way account for indirect correlations using a distance metric. I would not suggest trying to make your distance measure do all this work. Instead, I would use a simple but clear distance method, and then use a multidimensional analysis such as principal coordinate analysis (or similar) to understand the correlation structure of your compositional variables. This is a fairly standard workflow in metagenomics, for instance, which relies heavily on high-dimensional compositional data. I would recommend looking more into metagenomics statistical applications for inspiration if needed.
The choice does also depend to a certain extent on the specific application. In some fields certain methods are familiar and trusted even if there is not a strong empirical/theoretical rule.
You've asked a lot of questions about Pearson/Spearman correlations. Much has been written on the use of correlation analysis on compositional data, most of it cautioning against doing this kind of analysis without special treatment. This leads to difficulties in building meaningful association models for different variables. One method that I'm familiar with in metagenomics is sparse correlation, e.g. SPARCC. But there are newer methods that may be better, for example here and here.


text:
Upgrading guppy basecaller on Oxford Nanopore Gridion via PPA
I am trying to update the guppy_basecaller (current version installed is 3.2.10+aabd4ec)  on a ONT Gridion device. I am initially trying to update it via the nanopore PPA, as described here (requires login access to the nanopore community)
After I issued sudo apt update, I got a GPG error, saying that the public key was not found in the system.
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://europe.oxfordnanoportal.com/apt xenial-stable InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://mirror.oxfordnanoportal.com/apt xenial-stable-gridion InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch https://mirror.oxfordnanoportal.com/apt/dists/xenial-stable-gridion/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch http://mirror.oxfordnanoportal.com/apt/dists/xenial-stable/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Some index files failed to download. They have been ignored, or old ones used instead

I came across a possible solution shown here, so I tried to add the key to the system with the following command:
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys AD9D09AF2EBBA5A8

The terminal was busy for too long and, after some time I got an timeout error:
Executing: /tmp/tmp.EWnOg4cPJk/gpg.1.sh --keyserver
keyserver.ubuntu.com
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpg: keyserver timed out
gpg: keyserver receive failed: keyserver error

I then came across this post where one potential solution was to change the port, so I issued:
sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys AD9D09AF2EBBA5A8

But then I got:
Executing: /tmp/tmp.sfqFuJpVb2/gpg.1.sh --keyserver
hkp://keyserver.ubuntu.com:80
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpgkeys: key AD9D09AF2EBBA5A8 not found on keyserver
gpg: no valid OpenPGP data found.
gpg: Total number processed: 0
gpg: keyserver communications error: key not found
gpg: keyserver communications error: bad public key
gpg: keyserver receive failed: bad public key

Has anyone managed to update the guppy_basecaller via PPA, as described in the nanopore community page?
An alternative would be to install via the .tar.gz file, though I am first trying to upgrade the current version from the PPA to avoid conflicts with the current installation.

This is an ONT problem, not an ubuntu problem. Keys can expire over time, and need to be updated.
The keys for the Guppy basecaller are not present in the ubuntu repository; you need to get them from ONT. This is indicated by the installation instructions for Guppy (which you've already linked to):
sudo apt-get update
sudo apt-get install wget lsb-release
export PLATFORM=$(lsb_release -cs)
###
wget -O- https://mirror.oxfordnanoportal.com/apt/ont-repo.pub | sudo apt-key add -
###
echo "deb http://mirror.oxfordnanoportal.com/apt ${PLATFORM}-stable non-free" | sudo tee /etc/apt/sources.list.d/nanoporetech.sources.list
sudo apt-get update

Yes it turned out that the upgrade worked after repeating all these commands (which I thought it was not necessary because I had already added the PPA). There was, however, an issue with the release, so export PLATFORM=$(lsb_release -cs) fixed the issue


text:
Convert DNA raw data between a different formats
I have .vcf file with partly decoded genome provided by some service. How can I convert this file into other popular formats which may be used in other services? I've tried plink2 but unfortunately, rather I just can't use it correctly.
Also, where could I get a most completed list of a free DNA services?

Here's part of my DNA raw data with file headers:
##fileformat=VCFv4.2
##source=Genotek
##reference=hg19
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  vh6913
chr1    752721  rs3131972   A   G   .   .   .   GT  1/1
chr1    759036  rs114525117 G   A   .   .   .   GT  0/0
chr1    801536  rs79373928  T   G   .   .   .   GT  0/0
chr1    834830  rs116452738 G   A   .   .   .   GT  0/0
chr1    835092  rs72631887  T   G   .   .   .   GT  0/0
chr1    838555  rs4970383   C   A   .   .   .   GT  0/1
chr1    838665  rs28678693  T   C   .   .   .   GT  0/0
chr1    840753  rs4970382   T   C   .   .   .   GT  1/1
chr1    846808  rs4475691   C   T   .   .   .   GT  0/1
chr1    851390  rs72631889  G   T   .   .   .   GT  0/0
chr1    854250  rs7537756   A   G   .   .   .   GT  0/1
chr1    861808  rs13302982  A   G   .   .   .   GT  1/1
chr1    863130  rs376747791 A   G   .   .   .   GT  0/0
chr1    866893  rs2880024   T   C   .   .   .   GT  0/0
chr1    868404  rs13302914  C   T   .   .   .   GT  1/1
chr1    872952  rs76723341  C   T   .   .   .   GT  0/0
chr1    878331  rs148327885 C   T   .   .   .   GT  0/1
chr1    879911  rs143853699 G   A   .   .   .   GT  0/0
chr1    881627  rs2272757   G   A   .   .   .   GT  0/0
chr1    884767  rs67274836  G   A   .   .   .   GT  0/0
chr1    888659  rs3748597   T   C   .   .   .   GT  1/1
chr1    889238  rs3828049   G   A   .   .   .   GT  0/0
chr1    891277  rs77608078  C   T   .   .   .   GT  0/0
... and more than 636 thousand rows

Hi, you probably should be a bit more specific about which formats you need.

Yes, please [edit] your question and explain what kind of conversion you need. VCF file normally contain specific variants, so there are very limited options for "converting". I suspect you are thinking of a gVCF file with information for all positions of a genome, but we can't really help without more details.

@WouterDeCoster I need to get 23andme file format

@terdon I do not understand all these nuances about VCF, I just have a file with .vcf extension. Should I post a part of it for example?

Yes please. Post an example of your input file and the output you would want from that example. If you tried plink, please also tell us what you tried and how iot failed so we don't suggest the same solution.

@terdon thanks for your response, I updated question with data. Unfortunately, I don't know what is the format 23andme using, sorry. I thought this is kind of well known format. With plink I tried this solution but seems this not works this way for me, may be my plink version has differences (I use latest)

OK, that is a VCF file. But you still haven't told us what you want to convert it to. This is the format most commonly used to store variant information. What is it you need?

He said he needs to convert it to 23andme format. This is quite easy bcftools convert -c ID,CHROM,POS,AA -s SampleName -f 23andme-ref.fa --tsv2vcf 23andme.txt -Oz -o out.vcf.gz

@user438383 Your command converts from 23andme to VCF, not the other way around.

@RamRS yeah, I need vice versa. Could you help?

Isn't this all that's needed?  awk 'BEGIN{OFS="\t"} /^[^#]/{print $3, gsub("chr", "", $1), $2, $4$5}'  sample.vcf   > sample.23ame  (from a brief look at http://fileformats.archiveteam.org/wiki/23andMe#File_format )

The awk would work, sure, but one is better off using bcftools query -f

For completeness, that command would be bcftools query -f '%ID\t%CHROM\t%POS\t%REF%ALT\n' sample.vcf > sample.23andme but you might or might not have to strip of the "chr" part, depending on the further use I guess...


text:
Median string problem & multiple sequence alignment
I read about the median string problem as an introduction to the multiple sequence alignment, however none of the MSA algorithms used seems to be using the idea of finding the median string.
To my understanding from the consensus (median) string, it would be possible to retrieve the actual alignment by simply doing the pairwise alignment of the consensus string with each of the strings.
Why is this approach not used? Are there any algorithms that use it?
And related to this, what kind of tasks in bioinformatics/biology are median/consensus string algorithms used for?

Why is this approach not used? Are there any algorithms that use it?

As shown in this paper, the MEDIAN STRING problem is NP-Complete on alphabets of carnality greater than 3. While the problem isn't used directly in approaches due to the complexity, the idea of trying to create a consensus sequence between pairs of input and using that as building blocks is common. These are considered "progressive alignment methods".


text:
Are there any databases for gene co-expression or expression pattern clustering?
I am currently working on gene clustering based on co-expression pattern in mouse brain. The problem is I do not have some solid way to test my result. Are there any suggestions for databases indicating groups of genes that co-express (in mouse brain), or even better indicating some similarity metrics between each pair of genes based on expression pattern?
Also, I do want some databases for differentially-expressed genes ;)
PS: By "differentially-expressed" I am referring to a pair of genes whose expression patterns in a single individual differ.

How about:

http://www.informatics.jax.org/expression.shtml
http://www.emouseatlas.org/emage/home.php

Other potential databases of interest:

https://biokeanos.com/search?q=co-expression+%2Bmouse

Thanks for your great answer!  The first two seems to be developmental expressions, while the third database is very helpful!  This site looks quite good: https://gemma.msl.ubc.ca/


text:
Spanning deletion conversion
I have a vcf file created with GATK with spanning deletion https://gatk.broadinstitute.org/hc/en-us/articles/360035531912-Spanning-or-overlapping-deletions-allele- and I now need to compare this VCF with its gold standard using hap.py (https://github.com/Illumina/hap.py/blob/master/doc/happy.md#getting-started) but the spanning allele gives me the following error:
[W] too many AD fields at chr10:133017183 max_ad = 2 retrieved: 3 and I've read that is due to * in the ALT column.
So, is there a way to convert the spanning alleles into simple vcf files?
Thanks a lot,
Denise

Could you describe what problems you're having? You could try using rtg-tools vcfeval since it's designed to account for differences in indel representation between the truthset and callset.

Thanks a lot, I know I can use it but I would like to follow the best practises given by the Global Alliance for Genomics and Health Benchmarking Team (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6699627/)

Please [edit] your question and include some of the offending lines. Also make sure to include the relevant headers since the headers are where this sort of thing is defined. We can't really help without seeing the file.


text:
Create new range based from existing range (python)
I have a pyranges object that I am looking to iterate over. For each range, I would like to "resize" it to its centre (so start+end/2) and save that into a separate pyranges object. There will be additional processing associated based on what ranges already exist nearby, so I can't just use a resize function, I have to iterate over the object.
First, how would I go about manually creating a new pyranges object, then inserting one by one new ranges into it?
Second, when I try to iterate over my current pyranges object (using this line for myspot in pr.itergrs(hotspots, strand=False) :), it gives me this error: AttributeError: 'tuple' object has no attribute 'stranded'. This error occurs whether or not I have strands in my pyranges, and whether or not I specify strand=False. Why does this occur?

The data I am working with: a pyranges covering hotspots in the human genome (all chr1 for the moment)
+--------------+-----------+-----------+-----------+
| Chromosome   | Start     | End       | Score     |
| (category)   | (int32)   | (int32)   | (int64)   |
|--------------+-----------+-----------+-----------|
| chr1         | 154421820 | 154421821 | 26413     |
| chr1         | 30749469  | 30749470  | 14833     |
| chr1         | 55259503  | 55259504  | 9900      |
| chr1         | 14985097  | 14985098  | 7735      |
| ...          | ...       | ...       | ...       |
| chr1         | 247975751 | 247976497 | 1         |
| chr1         | 247977599 | 247981391 | 1         |
| chr1         | 248018012 | 248019454 | 1         |
| chr1         | 248177298 | 248178490 | 1         |
+--------------+-----------+-----------+-----------+

The hotspots are sorted by their score, and I'd like to go through and for each one, resize it to its centre, keeping the score, but if there is already a hotspot within a certain range (say +/-5kb), to discard it instead, so that I end up with clear chunks between hotspots. This is the "pseudocode" I was planning to do to implement this:
add centre of hotspot 1 to a new pyranges object
create ranges with score = -1 for 5kb flanking each side
for myhotspot in hotspots:
  resize to centre
  if not overlap with new pyranges object
    add centre
    create flanking regions
discard ranges with score<0

This may not be the most elegant solution, but I'll never have more than ~2000 ranges to process so it should work fine. Is this implementable with pyranges?

Quick and dirty (not to mention atrociously bad O(n)) solution:
import pyranges as pr
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=10000, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

def remove_worst_scores_until_no_overlap(gr):
    df = gr.df
    old_length = -1
    new_length = len(df)
    while old_length != new_length:
        df = df.drop(["Cluster", "Count"], axis=1, errors="ignore")
        df = pr.PyRanges(df).cluster(count=True).df.sort_values("Score", ascending=False)
        df = df[df.duplicated('Cluster', keep='last') | (df.Count == 1)]
        old_length = new_length
        new_length = len(df)

    return pr.PyRanges(df.drop("Cluster Count".split(), axis=1)).sort()

hotspots_no_overlaps = remove_worst_scores_until_no_overlap(gr)

Quick and dirty, but it works!

Answer to main question:
import pyranges as pr
assert int(pr.__version__.split(".")[2]) >= 93, "pip install pyranges==0.0.93"
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string("""Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
""")

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Could you please post an example of what data you have and what you want to end up with?

Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.

pr.itergrs is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.

To iterate over a pyranges you can do for k, df in gr: # ....
Anyways, to do what you want:
import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1

I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)

Code example pr.itergrs:
import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print("df:")
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -

Thanks for the pointers, this already helps a lot. I'll post my code in the morning (it's 7PM for me now)

Added the information you asked for :)

Thanks for the solution, do you think you could comment the function a bit more so I can be sure of what each part of it does? From what I can tell, you sort and cluster the df by score in f, then merge all ranges withing 5kb of each other, is that correct?

I collect all ranges within 5k and put them in separate clusters. Then the function f in apply merely picks out the row within each cluster with the highest score and discards the rest. I might have made a mistake in that if you have hotspots at 0, 4000, 8000 and 12000 then you would probably want hotspot 0 and 8000 while I discard all but one. Will think about the problem some more and update the answer tomorrow.

Yes I was wondering if it behaved like that. That's why I was thinking of creating "flanking" regions to better define which hotspots to discard, but if you find a more elegant solution, I'm all for it!

I've now added a method max_disjoint to PyRanges in 0.0.93. Thanks for making me realise there was a fundamental interval operation I was missing.

Surely "Dev adds new function to package because of you" should be a medal on bioinformatics SE, right? Thank you for your help in any case! I may be calling upon your services later depending on how the project advances

I do this for fun, but thanks. Feel free to accept the answer if everything worked as expected.

So the max_disjoint() works great but when I combine it with the centering you gave, it throws an error: ValueError: Buffer dtype mismatch, expected 'const int32_t' but got 'long long'. Do you want the full traceback?

Will fix. Could you post an issue on the pyranges github?

Thank you for all your help, but max_disjoint() has a behaviour that could be improved on I think: for now it seems to start at the first range on the chromosome, whereas it would be useful to be able to specify a column for it to sort on, like the score, before searching for overlaps.

Actually I'm not entirely sure what it uses to sort the pyranges on, seeing as even in your example, it's kept the 10001;10002;1 range as the centre of the second cluster

max_disjoint() is actually a well-known greedy algorithm: https://en.wikipedia.org/wiki/Maximum_disjoint_set#1-dimensional_intervals:_exact_polynomial_algorithm It starts by sorting on the end column. I will continue thinking about the problem


text:
How to identify specific motif for chip-seq data
For chip-seq data, I had analyzed the samples until motif discovery using both cluster profiler and Homer for each histone mark. I was wondering what is the efficient method to determine if each GO term or identified motifs has a biological meaning or not or if there is any way to filter out the genes/motifs for only specific cell types such as PBMCs, how to prioritize genes? I know random matches might even surpass the biologically relevant matches but right now I am not sure how to narrow down the motifs/genes? Your help would be greatly appreciated.

There is no easy answer to this without knowing what you are exactly doing and which question you want to answer. In general ChIP-seq is a good support experiment to answer a scientific question you have, deriving a question from it is hard. Please add some details, even though I am tempted to vote for closing as this is very unspecific and open-ended,  but lets see the details first.

To add more detail I have extracted PBMCs from infectious patients with Anthrax. I performed chip-seq for 3 histone marks: (H3K27ac, H3K4me3, H3K27me3).I'd like to know if the bacteria is affecting the epigenetic signature of PBMCs?
After peak calling and annotation for each marks, I have more than 2000 genes for each mark and more than 500 motifs called by Homer. I am trying to narrow the genes to those which I expect to see in PBMCs cells and those which have biological meaning. Right now I am struggling to find a way to filter out meaningful genes.


text:
SLURM script for running RSEM star fails
If you have worked on cluster with submitting jobs, please help me in resolving the following error.
I am running rsem pipeline to generate star index. In principal STAR requires more than 30G of RAM for human genome.
I have following script
#!/bin/bash
#SBATCH --mail-user=<your_email_address>
#SBATCH --mail-type=end,fail
#SBATCH --job-name="job01"
#SBATCH --mem-per-cpu=10G

rsem-prepare-reference --gtf ~/Genomes/Homo_sapiens.GRCh38.100.chr.gtf -p 5 --star --star-path ~/Tools/STAR/bin/Linux_x86_64/ ~/Genomes/Homo_sapiens.GRCh38.dna.toplevel.fa ~/Genomes/GRCh38_100

The process fails as soon as it starts to generate genome files using star. Below is the error complaining about the memory.
I have put 5 threads for star usage and 10G for each cpu.
here is the error
rsem-preref /storage/homefs/as21u994/Genomes/GRCh38_100.transcripts.fa 1 /storage/homefs/as21u994/Genomes/GRCh38_100
Refs.makeRefs finished!
Refs.saveRefs finished!
/storage/homefs/as21u994/Genomes/GRCh38_100.idx.fa is generated!
/storage/homefs/as21u994/Genomes/GRCh38_100.n2g.idx.fa is generated!

/storage/homefs/as21u994/Tools/STAR/bin/Linux_x86_64//STAR  --runThreadN 8  --runMode genomeGenerate  --genomeDir /storage/homefs/as21u994/Genomes  --genomeFastaFiles /storage/homefs/as21u994/Genomes/Homo_sapiens.GRCh38.dna.toplevel.fa  --sjdbGTFfile /storage/homefs/as21u994/Genomes/Homo_sapiens.GRCh38.100.chr.gtf  --sjdbOverhang 100  --outFileNamePrefix /storage/homefs/as21u994/Genomes/GRCh38_100
Feb 25 12:10:58 ..... started STAR run
Feb 25 12:10:58 ... starting to generate Genome files

EXITING because of FATAL PARAMETER ERROR: limitGenomeGenerateRAM=31000000000is too small for your genome
SOLUTION: please specify --limitGenomeGenerateRAM not less than 168632691637 and make that much RAM available

Feb 25 12:18:11 ...... FATAL ERROR, exiting
"/storage/homefs/as21u994/Tools/STAR/bin/Linux_x86_64//STAR  --runThreadN 8  --runMode genomeGenerate  --genomeDir /storage/homefs/as21u994/Genomes  --genomeFastaFiles /storage/homefs/as21u994/Genomes/Homo_sapiens.GRCh38.dna.toplevel.fa  --sjdbGTFfile /storage/homefs/as21u994/Genomes/Homo_sapiens.GRCh38.100.chr.gtf  --sjdbOverhang 100  --outFileNamePrefix /storage/homefs/as21u994/Genomes/GRCh38_100" failed! Plase check if you provide correct parameters/options for the pipeline!

Could anyone of you tell me how to fix this.
Thank you

Have you tried allocating more than 10GB of RAM? It seems like the error is related to memory..

I have tried till 50G. Don't know what is happening

Try to use the --mem= parameter to allocate a fixed amount of memory for the entire job, not mem per cpu. Your cluster probably has a default for that value which is probably too small for the job. For this either consult the cluster documentation or the admins.

Users have reported similar issues to STAR before: https://github.com/alexdobin/STAR/issues/129


text:
Snakemake fails because it cannot create logs inside a read-only Singularity Container
I am using Snakemake from inside a Singularity container. Because of that, the directory that Snakemake wants to log to cannot be created, and I get an OSError [Errno 30] Read-only file system: '.snakemake, more specifically the traceback
Traceback (most recent call last):
  File "/usr/local/bin/snakemake", line 8, in <module>
    sys.exit(main())
  File "/usr/local/lib/python3.8/dist-packages/snakemake/__init__.py", line 2562, in main
    success = snakemake(
  File "/usr/local/lib/python3.8/dist-packages/snakemake/__init__.py", line 525, in snakemake
    logger.setup_logfile()
  File "/usr/local/lib/python3.8/dist-packages/snakemake/logging.py", line 137, in setup_logfile
    os.makedirs(os.path.join(".snakemake", "log"), exist_ok=True)
  File "/usr/lib/python3.8/os.py", line 213, in makedirs
    makedirs(head, exist_ok=exist_ok)
  File "/usr/lib/python3.8/os.py", line 223, in makedirs
    mkdir(name, mode)
OSError: [Errno 30] Read-only file system: '.snakemake'

The Snakefile I am using is inside a Singularity-container (there's no way around that, I am on a hpc and don't have root, also using a --sandboxed container doesn't work), and I call it like this (I've also played around with binding other directories without success):
singularity run --bind $(realpath .) --bind ~/container.sif ~/container.sif test_simple --cores 4

The test_simple rule is a really simple one that only unzips a gzip, for testing - the problem however lies generally in Snakemake wanting to create the log at <snakemake_base>/.snakemake/log (see https://github.com/snakemake/snakemake/blob/bd5a9239b2684b57ab0143bb709e77afa1493cdc/snakemake/logging.py#L308), which in my case is inside the container.
It doesn't help to use an overlay --overlay overlay.img for some reason, and I've also tried using my own (empty) log-handle (--log-handler-script ~/loghandler.py) without any success.
Any ideas on how to circumvent this or disable the logging or write it to somewhere else?

I don't know what happens in your case, but I use Singularity containers with snakemake inside, and when it runs, the .snakemake folder is created inside the current working directory, not in a read-only folder. It might help seeing the container definition file.

Basic sanity check: do you run the container from a folder in which you have write access?

Can you write to the logs into /tmp/ ?

Or maybe don't run Snakemake in the container (e.g., conda install -c conda-forge -c bioconda snakemake doesn't need root).


text:
Dataset showing evolution of gene
As a way to exemplify how the hierarchical clustering of particular DNA sequences can be used to derive phylogenetic trees, I'm looking for a (even very small) dataset of the same genes (or other biological sequences) from different organisms (for example animals where the evolutionary tree is known). Do you have any pointers?
Thanks in advance

I don't understand what the problem is. If you just want a bunch of sequences of homologous genes from different organisms, it is very easy to get them from places like NCBI or Ensembl. What is giving you trouble here? Do you need help using the standard databases, or do you need help in selecting the genes of interest? Could you please [edit] your question and explain what you are stuck on?

If you're building a tree of different species, it's probably better to use protein sequences. A common resource is Uniprot.
If you still want DNA sequences you can turn to NCBI Nucleotide.
Some examples that give decent trees: p53, TOP1, HBA.
A common way to construct a phylogenetic tree is to start with a multiple sequence alignment (MSA).


text:
ATAC seq density calculation
This paper Fig2b

b Changes of chromatin accessibility in AMD (n = 14) relative to
normal (n = 11) in all retina samples. Each dot represents one
ATAC-Seq peak. Blue line in the left panel indicates average fold
changes of peaks with the same ATAC-Seq intensity. The percentage of
reduced peaks is shown under the density curve in the right panel.

How is the density calculated or what is being shown there? average ATAC seq signal or Fold change?

They are not referring to density of reads. It is the density of the fold change, the y variable in the plot to left of Fig2b. It tells you the distribution of the foldchange and 92% of the fold change is <0.
See this page for more details on density and histograms
Suppose your data is something like this (purely making this up):
set.seed(555)
fc = rnorm(1000,-0.5,2)
averATAC = -0.5*(fc - 0.5)^2 + 0.2*fc + 3  + rnorm(100,25,10)

You fit a density to it and you get the plot:
par(mfrow=c(1,2))
plot(averATAC,fc)
plot(dens$y,dens$x,xlab = "density",ylab="fc")

your answers are part of my phd work I mean those example which helps me to express my solution

how is the percentage being calculated from the density?

I think it's an overkill, but you can do sum(dens$y[dens$x<0])/sum(den$y)

"overkill" i will go for it

Just count the number of reads overlapping an ATAC-seq peak and divide it by the length of the peak gives you the read density. Each dot is a fold change, but to calculate fold change, average density from all the samples are taken, that is what the "average" is referring to.

my query is about the 92.3 % density what does it represent ? The fold change or the average ?

I see, it's the proportion of fold change falls below 0 on the y axis of the plot on the left

so is that the negative fold change they plotted?

it is a density plot of fold changes, not the density of the reads.


text:
How to identify genes from a genome assembly of C. Elegans?
I have two full genome assemblies for C. Elegans samples collected from two different geographical areas that I found on WormBase. These are in fasta format. I want to go gene-by-gene and compare the nucleotide sequences corresponding to the gene between the two samples. My goal with this is to count how many single-base differences there are in each gene. However, I'm not sure how to match the sequences in the assembly to known reference sequences (by the way is there a good resource to download all of these at once?) for C. Elegans genes.
My thought was to run BLAST between the reference sequence and each genome and find the match with the best score, if any. Is this a reasonable approach, or is there some better way to do this?

You don't want BLAST here, or at least not regular BLASTn. A better method would be to use tBLASTn to map the known C. Elegans proteins to the translated genome. However, I would recommend using a slightly more sophisticated approach:

Collect fasta sequences of all known C. Elegans proteins (amino acid sequences, not nucleotide). This should be easy enough to do on WormBase. For example, the protein sequences for C. elegans strain VC2010, assembly PRJEB28388, can be found here: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/species/c_elegans/PRJEB28388/c_elegans.PRJEB28388.WS279.protein.fa.gz

Use a tool that can model introns and exons to map a protein into a genome. I haven't worked in this field in more than 10 years now, so it is very likely that there are other programs around today, but back in the day, I would do this using either exonerate or genewise. With exonerate, the command would be:
exonerate -m protein2genome -n 1 -t assembly.fasta -q allproteins.pep > out

That will give you a list of best-hits for each of your input proteins and should accurately map the protein to the genome. Once you have that, you can extract the matches and start comparing.


text:
How can I get bedtools to tell me which genes are being expressed?
I'm trying to align the Acinetobacter baumanii genome to a genome. I've already done the alignment, and I want to use bedtools to see which genes are being expressed exactly. When I try running the following command,
bedtools coverage -a outCSEQSAligned.out.bam -b GCF_000746645.1_ASM74664v1_genomic.bed

unfortunately, it doesn't tell me which genes are being expressed at each place but instead tells me that the sequence aligns to the single Acinetobacter chromosome which isn't too helpful:
NZ_CP009257.1   282134  282254  GWNJ-1012:340:GW200702637th:1:1105:15519:10019/1        255     +       282134  282254  0,0,0   1       120,    0,4       120     120     1.0000000

I'm assuming that 282134-282254 are the positions on the chromosome that correspond to that alignment. I suppose I could extract these positions and figure out which genes they correspond to, but I'm wondering if there's an easier way to use bedtools to figure out which genes are being expressed. For reference, outCSEQSAligned.out.bam is my bamfile which I aligned the Acinetobacter genome to, and GCF_000746645.1_ASM74664v1_genomic.bed is a bedfile corresponding to the Acinetobacter genome.
Thanks for the help!

Did you do an RNASeq experiment?  It's not clear from what you wrote.

Bedtools is not the best tool for this. You can use featureCounts, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.
bedtools coverage is used to calculate the overlap between each read in BAM and each line in BED. It is not very straightforward to quantify gene expression. Plus bedtools does not account for exon in each gene and summarize them to gene level, like featureCounts does.


text:
Quantify low complexity regions in DNA sequences
I have fasta files with multiple sequences. They are reads that mapped to the genome outside of  probe-targeted regions.
From a quick perusal, they appear to be repetitive and have low complexity. Is there a way I can quantify their complexity and use that as an explanation as to why those regions were mapped to?

"Complexity" is also known as entropy. A quick and dirty way is to use a file compression utility - zip will get good compression for a repetitive sequence.
You can also look at the k-mer distribution for some suitable k: in a repetitive sequence it will be biased. Here is a tool that give a local plot: https://rdrr.io/github/vsbuffalo/qrqc/man/kmerEntropyPlot.html

I looked into compressibility as a measure of complexity. I found this article very interesting. They do not recommend it for DNA sequences.
https://sci-hub.do/10.1109/iwobi.2014.6913941

Use Picard EstimateLibraryComplexity.
From the docs:

Estimates the numbers of unique molecules in a sequencing library.

java -jar picard.jar EstimateLibraryComplexity \
     I=input.bam \
     O=est_lib_complex_metrics.txt

Although I mentioned I had fasta files, I was able to take the original bam files and pipe only those regions into picard with ´samtools view -b -L regions.bed sample.bed | picard....´

How long are your reads, and how many?  Some colleagues of mine wrote https://github.com/eclarke/komplexity to quickly score complexity in large numbers of short reads.  The score is derived from the number of unique k-mers normalized by the sequence length, so it's simple and fast, though the score inherently becomes less informative for longer and longer sequences (so, best for short reads).  We've used it for filtering repetitive or other low-complexity sequences in metagenomics projects.

This looks promising!
I don't want to run it on the unmapped per se as you can see in my question. But since all of the non-targeted regions are less than 300 bps, I think komplixity will be helpful.  I'll let  you know how it goes, and possibly award this as the accepted answer.

do you think you could ask your colleagues to take a look at their issues on github? I'm not sure it's being maintained.

I see what you mean, those issues aren't getting any attention.  I think the short answer for most of those is that it'd need new features to address them, like handling paired-end reads automatically.  (But for that use case you can at least filter forward and reverse separately, and take the intersection of what remains.)  If I can carve out the time I'll send eclarke some PRs to help with those though.

Thanks Jesse. That would be helpful.


text:
genotyping or variant calling in R
Are there any good options for calling variants or even just a pileup from R? R is not ideal for this, but I'd like to integrate with other functions.
I found deepSNV::bam2R which roughly does what I am looking for. Are there any alternatives?

You can always call out to the shell from R. Do you need the solution to be pure R, or could it be a library hosted on CRAN, bioconductor, or even conda? How are you controlling your dependencies?

@winni2k I was just looking for an R package. I am assuming the package is taking care of any dependencies.

Samtools/bcftools/tabix is available in R via Rsamtools.
SEQprocess also provides wrappers for different tools such as GATK and bowtie.
More utility packages are noted here.


text:
KeyError when getting features from a genbank file with biopython with some accessions but not others
I'm very new to python but i've been using it to extract the sequence of a gene from a genbank file. The issue is is that sometimes i'll get the output i want (prints the sequence to a file) and sometimes it will return a key error. This depends on which accession i'm using. Does anyone know why it sometimes might give a key error? I thought it might be something to do with the genbank records themselves, but they look pretty similar and the gene is there (in the gene feature qualifier). EG works with HG738867.1 but not AP019703.1. Here's my code -
from Bio import Entrez, SeqIO

gi_genome = 'accession'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
element = 0
for feature in record.features:
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
        report = 'record.features[%s]' % str(element)
        gene_sequence = feature.extract(record.seq)
        with open('output.fasta', 'a') as f:
            print('>' + gi_genome + ' ' + name, file=f)
            print(gene_sequence, file=f)
        break
    else:
        element = element + 1

Here's the traceback -
Traceback (most recent call last):
  File "/home/ubuntu/Documents/Git_Branches/Project_planning/Learning/In_progress/utils/data.py", line 11, in <module>
    if feature.type == 'CDS' and name in feature.qualifiers["gene"]:
KeyError: 'gene'

Process finished with exit code 1

Thanks in advance!

Good question.
feature.qualifiers is a dict, if the dict doesn't have that key it will throw a KeyError. The way your code works it expects feature.qualifiers to have 'gene' as a key every time feature.type == 'CDS'. However, it is possible to have feature.type == 'CDS' without a gene entry and in the case of AP019703.1 this comes up before you reach the gene acrA.
As simeple fix is to simply use get instead:
feature.qualifiers.get("gene",[])

Now returns a empty list if there is not key "gene", which also satisfies the 'name in' part of your if statement.
If you get stuck on stuff like this use more print statements, break up your if statements and look at the whole file (example below) to understand the structure.
gi_genome = 'AP019703.1'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
with open('AP019703.gb' ,'w') as fout:
    SeqIO.write(record, fout, 'genbank')

That's great, thanks!

This was amazing, Thank you!


text:
Create a NHX phylogeny format with r or python
I need help to add information in my phylogenetic tree on each branch.
I have two pieces of information at my disposal:

A phylogenetic tree in newick format :

(mammal:0.14,turtle:0.02,(rayfinfish:0.25,(frog:0.01,salamander:0.01)50:0.12)95:0.09);

A matrix with the edge number and information in the Information column:

Edge    Inf
7        0
6        0
5        2
4        1
3        2
2        2
1        2

Does someone have an idea in order to add these informations by creating a NHX format in python (with ete3 for instance) or with an R package ?
Here is the tree :

Brief answer Firstly this is not a final solution. The simple answer is to use TreeNode.add_feature() in ETE3 which would enable NHX format to be inserted into the tree. There are two possible ways to implement the solution, but I'm not aware of a 'reader'.  I'll return to this later under the section 'Question (again)'.
Note  I assume that the OP made a typo with 'witg ETE3' and probably means 'with ETE3'?
Interest However as an approach to molecular phylogeny it is interesting theory and has clear application. Not in a simple tree, but when the volume of data overwhelms optical analysis. I am not aware of this being deployed in phylogeny - not done a literature search though, but would readily automate processing of complex trees.
Theory The OP uses the term 'edge' rather than branch, so this is referring to depth first search (DFS) type style calculations. In Phylogenetics there are 'branches' and 'nodes' and these are distinct, in this line of theory, I think, 'edges' can be nodes ... possibly, it is possible an 'edge' connects two nodes, in which an 'edge' is synonymous with a 'branch'.
The difficult I have with the information is that edge 0 is not defined. In Phylogenetics edge 0 would probably be the root. It is possibly the 'missing branch' in the trifurcation at the base of the tree. In this line of theory however 0 may simply be a null category, because an edge is not necessarily the same thing as a branch.
Data structure The information present is paired. thus ... salamander is 7,0 ... frog is 6,0 ... turtle is 4, 2... We can work this out because the data structure is 7,6,5,4,3,2,1 for odd numbers.
7->0, 6->0, 5->2 and so on ...
Moreover, 7 is at the top of the tree whilst, 1 is at the bottom and the 'edge' is an arbitrary designation. Phylogenetics wouldn't label nodes in this format, labelling would usually be with respect to the base of the tree (root). If I look at this from my world-view there's no 0-> ... but 5 the ancestor of 7 and 6 goes to turtles.
Thus far the information isn't the DFS, its presented as the input, for example to calculate DFS, and can be used calculate all sorts of things in this style of graph theory.
Question (again) Back to the question, I am not aware of DFS style theory being used in trees, therefore I would not be aware of an interpreter that would read it. If the OP was just wanting to stash the data into a tree via  NHX it could be stacked onto the root lineage. There isn't a root lineage in this tree so I've just added a root at '1:0.01', i.e. a nominal distance to provide the OP with a place to stash the data there. This will read as a tree and the additional 1:0.01 doesn't affect the integrity of the tree - especially in this given phylogeny where the location of the root is not in dispute.
Idea 1
As described above the Edge data gets stored in the 1:0.01 inserted as an additional bracket at the end of the tree.
((mammal:0.14,turtle:0.02,(rayfinfish:0.25,(frog:0.01,salamander:0.01)50:0.12)95:0.09)1:0.01);

Idea 2
This is pairwise data and NXH format could insert the pairwise data for each edge (via TreeNode.add_feature()) for each 'edge' thus salamander would have 7-0 or 7,0 inserted using NHX. That is absolutely technically where the data belongs, but how this read later is more challenging.
If I get chance I'll update the question later with specific NHX output and possibly code to achieve it.

OP upvoted. Its a good question with interesting application.

The treedata book has a section on tree annotations that should work for you.


text:
Oxford nanopore promethion; does guppy3 separate fast5 into pass and fail? If not what does?
just a quick question about fast5 pass fail folders. I have been working under the assumption that the fast5 pass fail folders I was given as raw data from our ONT vendor came from them doing base calling with guppy 3.something. Is that right? (The person that worked for the vendor has since moved on so I can't ask them - so my question is really does guppy3 divide fast5s into pass/fail?)
I further assumed that given I am now doing my own basecalling with guppy4.4 that I should use all the fast5s (both pass and fail) as some of the fails might become pass with the improved NN. Again, am I right? If not could someone please explain?
I'm going to also put this question on the ONT community site, FYI
Thanks,
Liam

Only if they've enabled the --qscore_filtering flag when using Guppy; this flag isn't enabled by default when running Guppy from the command line.
It's a good idea to use all fast5 reads for recalling, because accuracy improves over time, so some that previously failed may slide into the pass folder.
FYI, I'd be much happier if people asked their nanopore sequencing questions here. It's a publicly-accessible website which is indexed by search engines, so makes it easier to find answers, and easier to see the answers when found (i.e. no login needed).
Note that when recalling, you set the input and output directories (I think via the -i and -s flags). These should be different directories, and ideally different output directories for each recalling attempt, to avoid duplicate file entries. The original fast5 files shouldn't be modified by recalling as long as the input and output directories are different.

Thanks. The issue I face is that I have now got fast5s of the same name in both pass and fail folders. I want to keep splitting fast5s on qscore. I'm concerned then when I do basecalling I will over write fast5s output.


text:
How to import a phylogenetic tree (in nexus format) with node support values in R
I have a phylogenetic tree in nexus format that was generated with RAxML. In the tree file, there are node support values. I want to plot this tree in R with the support values showing. However, when I import the tree file in R with the read.nexus function in ape, the node values do not appear in $node.labels. Can someone recommend another option in R?

RAxML, or raxml-ng? The former has both a bipartitionsBranchLabels file with a weird support formatting in the newick, where the support values are wrapped in square brackets ([ ]) and a bipartitions file, where they are normally formatted.

See page 49 of the RAxML v8.2.X manual

I recommend you this tutorial about ggtree. You will find the solution.
[https://guangchuangyu.github.io/ggtree-book/short-introduction-to-r.html][ggtree]


text:
Difference in Calculation through GEO2R and difference in Article
Here is microarray analysis where author uses GSE50632 series on an  Agilent microarray GPL17660 platform with 4 eligible samples which I guess are GSM1224991, GSM1224992, GSM1224993 and GSM1224994 among which two are healthy and two are patient samples.
However I found the difference in the calculations-
Through GEO2R analysis table it is-
ID  adj.P.Val   P.Value t   B   logFC   miRNA_ID_LIST   SPOT_ID
148000  1.86e-02    5.82e-04    3.47    -0.402958   3.91    hsa-miR-3195

However in the paper the calculated values are far different specially logFC value.

What could be the reason for this?

Thanks for answering your own question @VishalKumarSahu, could you please consider "accepting" your answer. It definitely helps with site stats

Here are two example of calculations of logFC validated in GEO2R analysis which depends if values are already log2 transformed. During GEOQuery analysis there is calculation to transform where expression set is checked for value distribution and log2 transformation is applied-
ex <- exprs(gset)
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
          (qx[6]-qx[1] > 50 && qx[2] > 0)
if (LogC) { ex[which(ex <= 0)] <- NaN
  exprs(gset) <- log2(ex) }

Example 1:
=> Sample Values are already log2 Transformed
n(TumorSamples) = 45
n(ControlSamples) = 46

a = average(TumorSamples) = 9.199053111
b = average(ControlSamples) = 4.988139348

logFC = |a - b| = 4.210914

Example 2:
=> Sample Values are not log2 transformed
n(TumorSamples) = 2
n(ControlSamples) = 2

a = average(log2(TumorSample1) + log2(TumorSample2) ... ) = (log2(84) + log2(21.8)) / 2 = 5.419287

b = average(log2(ControlSample1) + log2(ControlSample2) ... ) = (log2(624) + log2(665)) / 2 = 9.331306

logFC = |a - b| = 3.91202


text:
Find species from FASTA files
In a school project we got the task of finding the species from two FASTA files. We got some hints of what to look for, like it is a unicellular eukaryote and that we can use blast. We don't know if the two files are the same species.
So I wonder if there is a standard way of finding species from FASTA files using the terminal?
I have a background in informatics, so I'm not well known with bioinformatics tools.

Are you not allowed to use the web BLAST tool? That's what I use if I need to quickly find out the likely origin of a DNA sequence.
There are command-line ways to do the same thing, but I only use them if I need to search hundreds or thousands of sequences at once.


text:
How is the t-statistic value calculated in GEO2R or Bioconductor?
I was trying to calculate t-statistics using Python's scipy and numpy as np-
scipy.stats.ttest_ind_from_stats(
    NP.mean(x[['GSM1224991', 'GSM1224992']]),
    NP.std(x[['GSM1224991', 'GSM1224992']]),
    len(['GSM1224991', 'GSM1224992']),
    NP.mean(x[['GSM1224993', 'GSM1224994']]),
    NP.std(x[['GSM1224993', 'GSM1224994']]),
    len(['GSM1224993', 'GSM1224994']),
)

For example in the dataset given below, GSM1224991 & GSM1224992 are in first group and GSM1224993 & GSM1224994 form second group.
ID_REF,GSM1224991,GSM1224992,GSM1224993,GSM1224994
148000,6.436150368369428,6.499787040655854,4.430816798843313,3.0819099697950434

Here is the calculation given in GEO2R analysis for the same reference but while I calculate t test it is 0.029630766122806237 instead of -0.33768 something as shown below.
ID  adj.P.Val   P.Value t   B   logFC
148000  1   0.001426    -9.44   -0.33768    -2.513

How to calculate t-stat here? Or considering single row is wrong approach?

That is not trivial. GEO2R uses the Bayesian linear model-based framework limma https://bioconductor.org/packages/release/bioc/html/limma.html for the analysis. I doubt you can easily and meaningfully implement the relevant code in Python without a major effort and notable background knowledge. The user guide https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf has background information plus links to the relevant papers. Don't try to analyze things by hand, use expert software. If it is not significant, well, then it is not significant. A 2vs2 comparison is strongly underpowered, therefore lack of DEGs are not unexpected, especially for microarrays.

I am learning R and I will make transition to R for sure but I want to understand each and every step and function to satisfying my quest.


text:
What does genotyping data and beta values for SNPs mean for GWAS study?
I have a list of Genes with their alleles (RS IDs) associated with a specific condition. I have to find their genotyping data and the beta values from published studies. I looked through several studies and none of them had anything mentioned about the said terms?
Apparently its for calculating the polygenic score of the listed SNPs. Does that mean the genotyping data is the allelic possibilities and beta score as effect size? Effect size is the population size (n)?

The beta statistic is a common output field produced by programs that carry out genome-wide associations. There's a good explanation here.
I interpret a beta statistic as the degree to which a change in genotype from the reference allele alters the phenotype of interest. The association (and beta statistic) can be negative or positive. It should always be tied to a particular variant, usually called the "risk allele".
The beta statistic is related to the number of individuals found in each group with a particular genotype. In that sense it could be argued to be related to the population size (at least, assuming pre-grouping of populations by genotype), but that's not a typical interpretation.
The GWAS catalogue is a good place to start off looking for data associated with GWAS:
https://www.ebi.ac.uk/gwas/home
However, I think that's for summary statistics. If someone has asked you for raw genotype data, that is trickier to get, and tends to depend on the condition. Genotype data is personally-identifying information, and there are lots of issues associated with a blanket public release of such information, which is why it's so hard to find.
Without more information on the particular phenotype, I can't give much guidance, but useful places to look for most traits are dbGaP and UK Biobank. Both of those data sources need forms filled out and approved (to demonstrate an ethically-responsible research aim) in order to access genotype-level data.
There are a few public datasets around that include data uploaded by individuals (who have explicitly consented to their genotype data being public), but it tends to be trait-specific.


text:
STAR aligner multiple fastq files
I’m using STAR to align fastq files from SMART-seq2. I have raw data folder containing sub-folders with samples names the sub-folders each contain fastq file. How can I make a bash command in order for STAR to go in each of sub-folders and make alignment for the fastq file and retrieve the bam files in one folder?
Folder structure:
/plates/plate_1/rawdata/A1/A1_R1.fastq
/plates/plate_1/rawdata/B1/B1_R1.fastq
...
/plates/plate_1/rawdata/O24/O24_R1.fastq
/plates/plate_1/rawdata/P24/P24_R1.fastq

It’s 384 sample folders A1, A2, B3, F5,,, so on

Thanks I added the path structure

It seems the sequencing is single-end, and there is only one fastq per sample. If so:
staridx=/path/to/star/genome/index

mkdir /plates/star/
cd /plates/star/

files=$(find /plates/plate_1/rawdata/ -name "*_R1.fastq")

for i in ${files}
do

  sample=$(basename $i "_R1.fastq")

  STAR --genomeDir $staridx --runThreadN 8 --outFileNamePrefix ${sample}. \
    --readFilesIn $i --outSAMtype BAM

done

All fastq files will be mapped individually, and all bam files will be saved to /plates/star/.


text:
Brackets in .Rmd code blocks
I am making some R markdown documents for my stats class.  I've noticed I get an error in some functions when I have a code block with {r}  instead of just r.  Like so...
``` r
this code will be fine
```

``` {r}
this code will have an error
```

When should I be using the brackets for codeblocks in .Rmd and when should I leave them out?

Does it work without a space between the ``` and the {r}? It may be that the Rmd processor is getting confused about the spaces, but ignoring the bare r, and assuming the code is R by default.
The examples here don't have spaces:
https://rmarkdown.rstudio.com/lesson-3.html

Thank you I believe you may be right.  I am having other problems with my code.


text:
Build tree to calculate beta diversity
I am trying to calculate beta diversity* using the UniFrac distance on my data, each sample having OTUs values.
To calculate it, I saw it required building a tree with percentages that represented the OTUs.
However, I was unable to figure what percentage fit each OTU to build the tree from the instructions I found on the internet.
For example, I saw in this link: http://readiab.org/book/0.1.1/3/1#4.1.2,  that they build this tree:

Could you please explain to me how to build the tree?

This is the diversity between local sites and enables the total in a landscape to be calculated (gamma). The tree would therefore appear (presumably) to represent the diversity between sites.

The key thing we need to understand is the nature of the input. If this is Whitikar's beta this likely a species count. It could be something super clever like a 16S metagenomics between sites. I wouldn't have trouble making a tree (its UPGMA) but I'd need to know what the data represented. A, B, C are likely different foci

Can you give a little more information about what tools and data sets you are using? How did you generate the OTU table etc?
If the underlying data is nucleotide data and this is a phylogentic tree. It may well be a population genetic tree based for example in allele frequencies.
For classical phylogenetics then FastTree is a very simple tool to generate a decent phylogeny. If you have the sequences of each OTU, you can align them in MegaX with muscle, or simply using muscle from the command line if you prefer. You can then generate the tree with this command (presuming you are on windows):
FastTree.exe -gtr -nt < aligned_sequences.fasta > otus_tree.tre

I hope that helps!


text:
Extract features from fasta sequences and train the classifier
I am new to the bioinformatics field. I have positive and negative protein sequences for acetylation PTM. Now, I want to train a classifier, say SVM. What will be the next step? How can I convert these sequences into usable features? Any information or links would help.
>P31327_55|1|testing
KAQTAHIVLEDGTKMKGYSFGHPSSVA
>P31327_57|1|testing
QTAHIVLEDGTKMKGYSFGHPSSVAGE
>P31327_119|1|testing
APDTTALDELGLSKYLESNGIKVSGLL
>P31327_157|1|testing
LATKSLGQWLQEEKVPAIYGVDTRMLT

It may be worth checking what other folk do for motif classification, for example see http://elm.eu.org/elms/candidates.html#MOD_acetylation which lists the different motifs (different enzymes recognise different motifs). Although ELM uses metadata such as protein localisation to filter the motifs for much better results.

Q to others. a ML features (measurable property) is a different thing than the sequence feature present in a GenBank file (i.e. an annotation on the sequence, e.g. acetylation sites). This the second time I have seen this tag usage, is this something worth discussing in Meta?

Just a recommendation: Be careful when making your data sets.
P31327_57, P31327_119, P31327_157 are all annotated as "By similarity" to the mouse entry in Swiss-Prot https://www.uniprot.org/uniprot/P31327#ptm_processing, which means there is no direct experimental evidence in the human entry. If you do have publications that report experimental evidence for these PTMs, please don't hesitate to share them with the UniProt team.
I also recommend reading this UniProt help page about using UniProt data in negative datasets if you haven't already done so: https://www.uniprot.org/help/negative_datasets


text:
Calling GEOquerry error while loading the library
Note: this question has also been asked on Biostars
I started getting a strange error when trying to load GEOquerry with library(GEOquerry).
Error: package or namespace load failed for ‘GEOquery’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]): namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

When I run > devtools::install_github('r-lib/rlang') I get ...
Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]) :
  namespace ‘rlang’ 0.4.9 is already loaded, but >= 0.4.10 is required

Has anyone dealt with this before?  Restarting my R session has no effect.

rlang 0.4.10 is the current version at CRAN: https://cran.r-project.org/web/packages/rlang/index.html
Try to install it via install.packages, not via GitHub directly.


text:
Print specific columns in a matrix on the basis of sample id's in the header
I have a matrix file (expression.txt) in which the first column is a gene_id and from the second column the sample id's start. This matrix has 20,000 columns with sample id's corresponding to different tissue. I want to extract the samples of one tissue. Example:  SM-5GZZ7 SM-5GIEN SM-5EGGH are samples from brain tissue. So I want to separate the list of samples [(file=brain.txt;n=5000] corresponding to brain tissue from this matrix (with bash or R)
Name    SM-5GZZ7    SM-5EGHI    SM-5EGHJ    SM-5N9CS    SM-5GIEN    SM-5EGHH    SM-5EGGH
ENSG00000223972.5   0   0   0   0   0   0   0
ENSG00000227232.5   187 109 143 251 113 139 199
ENSG00000278267.1   0   0   1   0   0   0   0
ENSG00000243485.5   1   0   0   1   0   0   0
ENSG00000237613.2   0   0   0   0   0   0   0
ENSG00000268020.3   0   1   0   0   0   1   1

The file brain.txt looks like (Number of lines in this files = 5000)
SM-5GZZ7
SM-5GIEN
SM-5EGGH
SM-5GZXC
SM-5GZXO
SM-5GZY1
SM-5EGHK
SM-5EGHL

This should answer the question https://stackoverflow.com/questions/11098189/extract-columns-from-a-file-based-on-header-selected-from-another-file

Thank you soo much @user438383

Hi @user438383. I tried the script by dogbane, it didn't work

Example in R. If you have a file of column names, names.txt, such as:
mpg
cyl
disp

and a data.frame like mtcars
> mtcars
                   mpg cyl  disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
 [ reached 'max' / getOption("max.print") -- omitted 23 rows ]

Then you can select the columns by:
library(dplyr)
list_names = read.table("names.txt", header=F) %>% dplyr::select(V1)
> mtcars %>% dplyr::select(list_names$V1)
                     mpg cyl  disp
Mazda RX4           21.0   6 160.0
Mazda RX4 Wag       21.0   6 160.0
Datsun 710          22.8   4 108.0
Hornet 4 Drive      21.4   6 258.0
Hornet Sportabout   18.7   8 360.0
Valiant             18.1   6 225.0
Duster 360          14.3   8 360.0
Merc 240D           24.4   4 146.7

Thank you @user438383 for the help

There's a million different ways to do simple taks like this. Pick one and learn it:) Python is my go to as it can do pretty much everything you'd ever want. I saved your example as csv with no spaces called expression.txt.
import pandas as pd
df = pd.read_csv('expression.txt')
df[['SM-5GZZ7','SM-5GIEN','SM-5EGGH']]

SM-5GZZ7    SM-5GIEN    SM-5EGGH
0   0   0.0 0.0
1   187 109 143 251 113 139 199 NaN NaN
2   0   0.0 0.0
3   1   0.0 0.0
4   0   0.0 0.0
5   0   0.0 1.0

I too prefer Python :)

Thank you for your answer. But I prefer to use bash. I have used awk in past to do this task (but only too get 10-12 columns). In this case, I have to extract 5000 columns based on a list in another file named (brain.txt).

@Liam McIntyre: Can you something in bash or R? and I have to extract nearly 5000 columns from the mastermatrix.txt file. All the columns to be extracted are in a file brain.txt.

I don't know how to off the top of my head sorry. I see that you've been given solutions in both bash and r. Now it's up to you to do some work to make them fit Ur problem.


text:
Error during File conversion
Conversion of a multi-entry sdf to a multi-entry pdb file with open babel results in the loss of the compound names in the resulting PDB file. Why do you think this might occur and how can I resolve this issue.

Add more information so people can help you better. Are you using the GUI or the command line version of ob? If GUI, what are your steps, and if CLI, what are your commands?

@RamRS I am using OB GUI. I opted for sdf as the input format for file conversion and pdb as the output file format for the conversion.

http://www.cheminfo.org/Chemistry/Cheminformatics/FormatConverter/index.html just to be sure people can try it

ouch multi molecules doesnt work online