text:
local BLAST error: BLAST Database error: Error: Not a valid version 4 database
I want to use a local blastn command to BLAST a multifasta file ~700 sequences using following command: blastn -db nt -query fasta_all.fasta -num_alignments 2 -out fasta_blasted.txt but I receive error :BLAST Database error: Error: Not a valid version 4 database.
I use a local nt database which I downloaded today (11th of May 2020) by running update_blastdb.pl --decompress nt [*]. I downloaded this but I did not do any post-processing and just the "raw" download files are in the folder (see image below which shows only a part of the files)

The blast+ I downloaded and installed is ncbi-blast-2.10.0+-x64-linux.tar.gz from https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/
So, both database and blastn version should be the latest. It probably has something to do with either the version of the database or the software, but as it both are the latest version, I find it strange that these would not be compatible. Any ideas what causes this error?

first I recommend you to check your downloaded data with md5 to be sure there was no problem with that step. How many files did you download?
Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called is a version corresponding to the database you are using.

Thank you for your reply. I downloaded 24 files in total, from nt.00 till nt.23. For every nt file there are following extensions available: .nhd .nhi .nhr .nin .nnd .nni .nog .nsq. For every nt file there is a nt.xx.tar.gz.md5 file available. Does that make sense?

Try to run md5sum nt.xx.tar.gz.md5 for each xx in 00 to 23. Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called from your latest version folder.

I think this last issue you adres could be the problem. When I run blastn -version it returns Package: blast 2.6.0, build Jan 15 2017 17:12:27 so not 2.10 which I downloaded. Once I also downloaded blast with conda so perhaps that is the problem.

Thomas, indeed this was the problem, very good suggestion. I removed my BLAST through conda remove blast to remove the old version. I followed the exact steps of https://www.ncbi.nlm.nih.gov/books/NBK52640/ but this did not lead to a correct installation. The command blastn was not found even after right configuration (check by echo $PATH). Now I recall that I install blast also through conda because the description of NCBI did not work. I now managed to install blast 2.10 by following these steps: https://ngsprimer.wordpress.com/. I receive no error now. Thanks!

Happy to hear it solved the problem! good luck!

I run into this problem and indeed it was mismatch of database and blast versions. So other don't have to fish it out of the comment section to, I added it to the answer (you as much your words as possible). Feel free to adjust it, but that's the solution that should be accepted I think.


text:
How do I generate a list of post-synaptic gene markers that will guide me in my search for these markers in intestinal stem cells?
I would like to generate a list of gene markers that represent post-synaptic related genes. I will then use this list to search through single cell sequencing data, in order to assess whether certain cell types express post-synaptic related genes. I think there might be a way to get this list of gene markers from NCBI.
Thank you!

something like this ? https://www.gsea-msigdb.org/gsea/msigdb/geneset_page.jsp?geneSetName=REACTOME_NEUROTRANSMITTER_RECEPTOR_BINDING_AND_DOWNSTREAM_TRANSMISSION_IN_THE_POSTSYNAPTIC_CELL

Yes! i think so. Thanks for sharing. I didn't know about this database!

source: Stack Exchange

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

sort: votes


text:
Tools for creating gene coexpression networks
What are some good tools (in Python if possible) for creating gene co-expression networks from RNA seq data?

After looking around, it seems like the Camoco library is suitable.

It would be nice to hear from the R crowd

it depends on the size of your data. If you have enough RAM available, then you may directly compute a correlation-based distance matrix using the cor() function (with method="spearman" or method="pearson"). But you may also give a try to Mutual Information distance using the Parmigene (https://cran.r-project.org/web/packages/parmigene/index.html) R package or partial correlations with the Corpcor (https://cran.r-project.org/web/packages/corpcor/index.html) R package.
If you have a large expression matrix, remember to use the excellent data.table package and its fread() function.

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


text:
Separate multiple sequence into different file, one sequence per file
I have a file contain multiple sequence, and I want to separate them by "gene:" into different file.
example:
example.fa
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

How to get results like this? One sequence per file, and after > is gene id, not cds id
KQK21959.fa
>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

KQK14546.fa
>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

A simple Biopython solution:
from Bio import SeqIO

for record in SeqIO.parse("example.fa", "fasta"):
    SeqIO.write(record, record.id + ".fa", "fasta")

Note, if you have multiple sequences with the same id, only the last sequence will be kept in the final output files

Thank you! It works excellent.

There's an app for that! The exonerate package (see here) includes a nice fastaexplode command that does exactly what you're looking for and will be much, much faster than anything we give you here since it is written in C and not a scripting language. Just install exonerate (you can do that with sudo apt install exonerate on Debian-based systems) and then:
fastaexplode file.fa

This will create these two files (file.fa is the original input file from your question):
$ ls
file.fa  KQK14546.fa  KQK21959.fa

That really depends on the tools you are using. If anything will do you can use the following python snippet.
If you know that each sequence is on a single line and your file isnt gigantic use this
EDIT
Updated so the file isn't read to memory at once to address @terdon point. Although, the above solutions are probably much better.
with open( "seqs.fasta" ) as f:
    out = ""

    for line in f:
        if not line.strip(): continue
        if line.startswith( ">" ):
            if not type( out ) == str: out.close()
            gene = line[1:].strip().split()[0]
            out = open( f"{gene}.fa", "w" )
            out.write( line )
            continue
        out.write( line.strip() )

This can also be done using awk. This is harder to read, but probably faster than python:
awk 'sub(/^>/,""){
        FILE=$0".fa";
    print ">"$0 >FILE;
        next
    }{print > FILE}' file.fa.

This will define a new output file FILE at each header line, whose filename is the header minus the > character, suffixed with .fa. And all lines in the file are appended to the current output file.
EDIT: Changed to a faster version found on: https://crashcourse.housegordon.org/split-fasta-files.html#using-awk-instead-of-csplitsed

This worked for me today:
# Load required packages
BiocManager::install("Biostrings")

library(Biostrings)

# Read the multi-genome FASTA file
multi_genome_fasta <- readDNAStringSet("*.fasta")

# Extract individual genomes and save as separate FASTA files
for (i in 1:length(multi_genome_fasta)) {
    header <- names(multi_genome_fasta)[i]
    sequence <- as.character(multi_genome_fasta[i])

    # Combine header and sequence
    combined_data <- paste(">", header, sequence)

    # Save the combined data as a separate FASTA file
    file_name <- paste0(header, ".fasta")
    writeLines(combined_data, file_name)

    cat("Saved", header, "as", file_name, "\n")
}

cat("Separation completed.\n")

There are many ways to do it. Here is one way to do in Python.
import os

with open("input.fasta") as f:
    for line in f:
        if line.startswith('>'):
            record_name = line.lstrip(">")
            outfilename = f"{(record_name.rstrip())}.fasta"
            if not os.path.exists(outfilename):
                out = open(outfilename, 'w+')
            else:
                print('File already exists!')
        out.write(line)


text:
Is it important to filter out poor quality cells before performing an integration analysis on single cell RNA sequencing data?
In order to perform an integration analysis of single cell RNA seq data, is it important to check the percentage expression of mitochondrial genes of cells as well as the feature counts to exclude dying or dead cells as well as multiplets from further analyses?
I ask because the standard Seurat work flow to perform integration as outlined in Integration and Label Transfer does not mention filtering steps prior to integration.
In this publication, Eleanor Alenzi et al. 2018 also did a canonical correlation analysis with Seurat but did not mention steps to filter out poor quality cells before doing the integrated analysis. This second publication does same. Both publications do not mention if their visualisations or differential gene expression analyses are done on integrated (discouraged) or raw uncorrected RNA data (recommended).
Following the styles of the publications linked above (many more do same), is it also efficient to perform integration without filtering out poor quality cells?

It is absolutely necessary to remove low quality cells: In the case of CCA (and this applies to other "integration" or "data alignment" methods as well), one would need to use "anchors", basically same type/state/kind of cells from the samples to be "integrated" and are used to "align" the different samples. In a scenario where you have a set of cells whose distinguished feature is high mitochondrial content due to apoptosis irrespective of their origin (epithelial, T cell, ...) and such cells are picked up as so called anchors, the whole integration effort will be off. Even when not using integration methods, bad quality cells quite often form their own clusters (highly skewed biology) not contributing anything to analyses and are often removed at the downstream steps anyway.
And when it comes to published work, unfortunately it comes down to the backgrounds of the reviewers. I tend to think that there are way too many single cell studies than experts with required bioinformatics and statistics so that a thorough review can be conducted. And lack of a paragraph on QC does not mean no QC has been performed anyway.

I usually run the usual standard Seurat - Guided Clustering Tutorial protocol to filter out poor quality cells (cells with high mitochondrial content and as well as multiplets). In the papers cell filtering as a quality control process is not mentioned at all. If I understood you well, does it mean that the presence of dead cells has no effect on the clustering process or clustering outcome?

Depending on what makes the cells bad quality and how different these bad quality cells from the rest, the overall clustering pattern might be effected (highly or partially) or not. I would nevertheless remove the bad cells and also would use SCTransform (or its scran counterpart) over crude log normalization in the downstream.


text:
How to understand and analyse RNA-seq data (for a beginner)?
I am trying to understand expression of a certain protein across Pseudomonas species. I downloaded an SRA file from NCBI and converted it to a fastaq file. I am not able to understand how to interpret and analyse this file. I want to get the mRNA sequence and see which genes are being expressed and find the mRNA of interest.
How do I proceed with getting the mRNA sequences and then for mapping them to genome?

For some studies, in GEO they also include the processed counts. You can also see whether it's the case for your study

SRA files are raw sequencing data that have been converted to the NCBI storage format. These can be converted to fastq files, e.g. with fastq-dump from the SRAtoolkit.
Beyond that I suggest you start with the Hitchhiker's GUide to RNA-seq and the DESeq2 workflow to get a background. Be sure to extensively browse the web for vignettes, guides and tutorials. YOu will find helpful resources to get started.


text:
Should you filter GWAS hits with high standard error?
I'm trying to figure out if I should be filtering out GWAS hits that have high standard error and I'm not quite sure what to do.  It seems like it might not matter, because the standard error is used to calculate the t-statistic, which is then used to calculate the p-value.  So in a way it's already built in.  But reporting SNPs that have very high standard error doesn't seem quite right.  What's the best way to handle this?

I am curious about the relationship between your p-value (or effect sizes) and standard error. I would expect the significant signals to have smaller stand error compared to the non-significant, background signals. If this is the case, there is no need to do filter.

That hasn't been my observation in looking at the raw UK Biobank GWAS results. There may be some correlation between effect size and standard error, but it's not a universal rule. I have seen results with large standard error, and also large effect sizes.

I agree. But if you look at genome-wide significant SNPs, are their std really large? Do you have rigorous filtering criteria before GWAS? I am also wondering what is your purpose of doing the GWAS? If it is just a standard GWAS or genomic correlation or LDSC I think it is fine without filtering? I do not remember seeing any paper doing that.

I'm not doing the GWAS, I'm only looking at other people's results. It's not typical for people to report SNPs with large SE in their tables, because that means the p-value will be quite high (and most people concentrate on p-value above all else).

I think the effect size is more of an issue than the standard error.
If the standard error suggests that the effect direction could change sign, then it might be a good idea to filter it out. Otherwise, if the effect size remains large (or at least positive) even after accounting for a large error, then [to me] that would attach more weight to the likelihood that the association is valid.
I did a very quick analysis of the results of one UK Biobank study on Twitter, where I used the statistic of "the least extreme value in a 95% confidence interval" in a Manhattan plot instead of the more commonly-used p-value:
https://twitter.com/gringene_bio/status/1207617723586371584


text:
Importing a WGCNA co-expression network into a networkX graph in Python
Is there a way to generate a WGCNA co-expression graph in R, and then to import that graph into Python such that I can represent it using a networkX object?

You can use exportNetworkToCytoscape or exportNetworkToVisANT to generate and egde list file with the two nodes at the ends and weight for each edge. Hopefully there's a way to import that into networkX.


text:
How to align Illumina DNA seq reads to CDS?
I have the CDS sequences for my genome. I want to align Illumina DNA seq reads from related fungal strains to these. This is to identify SNPs in just the coding sequences. I am not sure if I should use a DNA seq aligner like BWA and Bowtie, or to use a split aware aligner like STAR.
I have the full genome. But I wanted the SNPs for each of the CDS sequences, so I can create consensus sequences for these CDS for the related strains.
Can you provide examples or links to get alignment from CDS region or using VCFtools to get SNPs and INDELs from the CDS region? I have tried looking but haven't found anything useful. Thank you.

Answer from @fabio-marroni, converted from comment:
I suggest that you align to the full genome and after this, you only extract the SNPs in the CDS. This procedure is commonly used in all exome sequencing projects, so you can find a lot of instructions software. You can choose between extracting the regions from the alignment (using samtools), or keeping everything, calling the SNPs and then only select those in CDS using vcftools.
For looking at the CDS region, you could try bedtools, as suggested here.


text:
Gene Ontology equivalent - higher specificity?
I've recently discovered OrthoVenn2, which I'm using to compare proteomes and extract likely clade-specific genes. I'm enjoying software, and one of the things it offers is a breakdown of protein clusters in terms of Gene Ontology (example from a mock dataset). However, I find the Gene Ontology terms a bit too general, and would appreciate more specificity for my data.
I guess something more like the 'Subsystems' figure on MG-Rast datasets? In an ideal case, I'd enjoy something interactive where I could explore different levels of annotation (e.g. if I click on 'transferase' in the pie-chart, it would be nice to see 'methyl-trasferase', 'formyl-trasferase' etc), but I understand that might not exist and I'd be happy with something that requires a bit of command line usage too.
I'm working with a relatively medium number of genes (50-400) at one time.
Thank you for your time!

How are the GO terms assigned? If it's just showing you the terms that are overrepresented in your group of proteins, then those are probably the most specific terms that are overrepresented and any more specific terms would not be overrepresented and hence are not shown. GO itself can be extremely specific, but that doesn't mean your data will be.

@terdon Thank you. It looks like the software is using Swiss-Prot to cut down on computer memory, which might leave out some terms. I guess it might help if I ran the sequences I got there through another GO-finding software?

What? No, swiss-prot is just the well annotated version of trEMbl, it's a high confidence set of proteins and protein annotations. Nothing wrong with using swissprot and I don't see how that's relevant. My point is that if you are looking for overrepresented terms, then you will never find the most specific terms, you will find broader ones that are shared by many proteins. I don't think this is a limitation of the tool you are using, but of your data: why do you think that more specific terms exist in your case?

@terdon Because I ran the same sequences through BlastKOALA (KEGG-based) and the results were more in line with what I had in mind (broad terms like 'Metabolic pathways' and 'Lipid metabolism', but also the possibility to zoom in to subgroups (such as 'Glycerolipid metabolism' or 'Glycerophospholipid metabolism', and even individual enzymes). I guess I'd like something similar to that.


text:
When should I apply a multiple comparison adjustment?
I am performing a RQL analysis ("http://www.esapubs.org/archive/ecol/E095/002/suppl-1.pdf").
This analysis meassures the relationship between different traits of three different matrix showing if there is a positive or negative relationship. This RQL analysis make a fourth-corner method.
Well, this method have the option to  make a p-adjust method to solve multiple comparison problem. However, my variables are independent and I belive that I should not apply to my dataset. Am I right? Or always that two or more matrix are compared I should apply a multiple comparison adjustment?
Thanks in advance.

The p-value adjustments corrects the inflated number of significant results when performing many significance tests.
If you have a test with 5 % false positive rate and but you run it several times the chances of getting a false positive results in any of these tests is far greater than 5 %. For example repeating the test with 5 % false positive rate 20 times means that there is a 60 % chance of at least one false positive result.
You can see this in the tutorial you have linked to in your question, where the number of significant differences goes down from 26 without correction to 18 with correction.
Therefore, if you are performing multiple comparisons you should do a correction. The problem of inflated p-values will occur when you are performing many tests even if the tested variables are independent.


text:
VCF-merge fails due to tabix not producing .tbi files
I'm trying to use vcf-merge to combine 2 exome capture vcf files (~250K and ~330K in size) before trying it on all 96 samples. I'd appreciate any advice on the best way to do that! I've detailed what I've tried below. My issue seems to be with using tabix to convert the files to .tbi format.
Step 1: BGZIP
So far, I've zipped the files without issue:
bgzip sample1.vcf
bgzip sample2.vcf

which produces:
sample1.vcf.gz
sample2.vcf.gz

Step 2: TABIX
When I try to use this command:
tabix -h -p vcf sample1.vcf.gz

the stderr is:

Region 536999277..536999278 cannot be stored in a tbi index. Try using
a csi index with min_shift = 14, n_lvls >= 6. tbx_index_build failed:
sample1.vcf.gz

Using the -C option which works:
tabix -C -h -p vcf sample1.vcf.gz
tabix -C -h -p vcf sample2.vcf.gz

which produces:
sample1.vcf.gz.csi
sample2.vcf.gz.csi
Step 3: VCF-MERGE
When I use this command:
vcf-merge sample1.vcf.gz.csi sample2.vcf.gz.csi > out.vcf.gz

the merge fails, and I get this stderr:

Broken VCF header, no column names? at /usr/share/perl5/Vcf.pm line
172, <ANONIO> line 1. Vcf::throw(Vcf4_2=HASH(0x5645a760bc38),
"Broken VCF header, no column names?") called at
/usr/share/perl5/Vcf.pm line 867
VcfReader::_read_column_names(Vcf4_2=HASH(0x5645a760bc38)) called at
/usr/share/perl5/Vcf.pm line 602
VcfReader::parse_header(Vcf4_2=HASH(0x5645a760bc38)) called at
/usr/bin/vcf-merge line 183 main::init_cols(HASH(0x5645a761f438),
Vcf4_2=HASH(0x5645a760b248)) called at /usr/bin/vcf-merge line 279
main::merge_vcf_files(HASH(0x5645a761f438)) called at
/usr/bin/vcf-merge line 12

If I exclude the tabix files, the merge still fails and the stderr says:

The column names not tab-separated? Could not load .tbi index of
sample1.vcf.gz. The command exited with an error. Is the file tabix
indexed?

Any suggestions on how I can correctly merge these files?

Cross-posted on Biostars: https://www.biostars.org/p/437957/

Would it be possible to share the files with us? Do you really have a chromosome that's 536999278 bases long? What species are you working on?

Hi, I'm working with conifers which have very long chromosomes. As discusses below, the tbi cutoff is ~536,870,912 bp, so csi is needed.

You should go with the CSI index as it seems like the BAI indexes cannot handle chromosomes longer than 2^29.
Your problem is that you're trying to merge the index files using vcf-merge. You should merge the VCF files, not the .csi index files.

Do you have a reference for this limitation? Is it mentioned in the docs somewhere?

Yeah, I found it in a GitHub issue for the mosdepth/samtools repos: https://github.com/brentp/mosdepth/issues/41 https://github.com/samtools/samtools/issues/241

This has come up on the mailing lists and forums often enough over the years to be part of the common lore. It's spelt out in the BAI spec in SAMv1.pdf §5.1.1 and will soon be mentioned in the samtools index manual page (PR #1248).

Thanks for the advice, and sorry for the cross-post. You're right, CSI must be used if indexing is necessary since all of my chromosomes are larger than 2^29.


text:
Standalone Blast+; Automating searches, Blastn formatting, interpretation
I am struggling with using NCBI's standalone blast, particularly the blastn feature.
I have a large number of nucleotide sequences that I want to use blastn to ID. The goal is to retrieve information from the blast program and eventually import this information to excel. I want to know the capabilities of, or the easiest way, to format the output. I am running a remote search to NCBI's nucleotide database and want to filter and organize the information retrieved as much as possible. The ultimate goal is automate a way to run a remote search to NCBI's nucleotide database, and to work with the results in excel.
I am extremely new to this and struggling to find the best path. From what I can gather, many prefer using external programs, scripts, parsers, interpreters etc to organize their results. I am okay with doing this as well, but would still like to format the output of my blastn search as much as possible beforehand. It is returning a lot of data to me that I don't need.
Can I automate this in excel, and just use VBA/macros to further organize my data? It seems some people prefer to just import the output text file to a spreadsheet and work from there. I have also seen options like BioEdit, eutilities, APIs... I'm not sure which route I should invest my time into figuring out.
Thanks so much for lending your time and expertise!

You just have to change the output format using the -outfmt parameter. To have an output in a tabular format just use -outfmt 6. The resulting table have 12 columns containing the following informations queryacc subjectacc pident length mismatch gapopen querystart queryend subjectstart subjectend evalue bitscore. You may restrict the number of hits per query using -max_target_seqs 5 for example.

Hey Thomas, thanks for the reply.
I found some of these features previously but unfortunately it does not include nearly everything I want to include. The specifiers are missing some crucial information, and I would like to use more filtering options if possible. I want to somehow write something that will choose the best choice based on my criteria (probably a lot of if, else statements) to ultimately have 1 best choice for each query. I have over 30,000 queries, so I really want to minimize the amount of data that is returned to me.

I haven't found a way to have the query sequence returned in my output, nor the binomial of the subject. Since I will be having so much data returned to me, what is the best option I can use for choosing the best choice based on my criteria? Like setting thresholds and if/then statements. Also, unfortunately many of the default specifiers included for that format are not necessary for me, and many bits of info that I do want returned don't seem to have specifiers.. Its very confusing!

Hi @Sarah, I recommend you to keep the output format as a tabular separated file. You may use a -max_target_seqs 1 but it is not recommended for large searches ( https://academic.oup.com/bioinformatics/article-abstract/35/9/1613/5106166?redirectedFrom=fulltext). I think you should think about using R for example to further process your output table a write a loop to get the best hit for each and integrate if/else loops. Then you should use the  blastdbcmd to get the subject or query sequence.


text:
Blast results filters
I performed a blastn search of NGS data against ssRNA database download from Internet, with a expected value 10-4. The size of NGS data reads is of 125 bp.
I have analyzed the blast results of the libraries. These libraries are from different habitats. The most difficult habitat, oakwood, produced a blast results with a maximum length of 51 bp. These results have a qcov of 100 (in most of the cases).
I would like to know if these low length sizes comes from degradation?
Should do  I filter these samples with a length = 125 bp with the aim to obtain results with better quality or with a qcov of 100 is enought to get reliable results?
Thanks you in advance.

I think you should maybe rely on the hitscore, because it integrates both overlap length and % of identity. In addition, you may have a look at MMSeqs2 as an alternative to BLAST. It is much faster, and has tutorials to blast ssRNA db such as SILVA (https://github.com/soedinglab/mmseqs2/wiki).

What would "better results" mean? What is wrong with what you get now? If your queries have a maximum length of 125, why would you search for >= 125? How can you have a hit with a length greater than your query sequence? Please [edit] your question and clarify.

Could you please explain your RNA extraction protocol from field collection to lab?

The RNA extraction was extracted with a CTAB kit from tissue plant of different plants species. This plant species comes from the habitat wilder. My concern is that the blast results gave to me hits with a high qcov and pident and low length.

could please give more information? which ssRNA database are you using? and how your reads have been preprocessed (trimmed)?

Any chance these hits are viral species? This would provide a clear explanation about the phenomena.

yes, this hits must be a viral species, the database was a ssRNA database downloaded from NCBI. Moreover I trimmed the adapters with cutadapt. The problem is that oakwood habitat is the unique habitat that show this phenomena, the other habitats show longer length in the blast results.

I had a think about this one and am stumped. I thought there could be a connection with acute oak disease (AOD) or else chronic oak disease (COD), which are emergent infections threatening the population. However, they appear to be fungal and bacterial infections rather than viral.
Single stranded RNA plant viruses are very common, more so than in mammals, but to achieve the results you are obtaining the trees would need to be really sick if the predominent RNA was viral. Short viral RNAs are produced regularly by viruses as defective RNA and its a major area of research, considering a mechanism of defective interference,
In context it would be weird if this RNA swamped the tree/plant, you need an actively replicating virus for this stuff to appear. So I don't think its either viral or defective viral RNA.
The antiviral response in mammals can involve RNase and this could explain your degradation. Essentially, mammals (probably vertebrates) kick-back hard when virus appears. However, I couldn't really find evidence of this in plants.
I don't think the results you obtained have a clear biological explanation in pathogenesis. The best explanation is it was wet lab-technique, e.g. library construction resulted in overdigestion/fragmentation. The only thing I can suggest is template assembly using a normally distributed plant/oak RNA library.
I was probably overthinking in forging a link with an emergent (tree) pathogen.

!

Thanks you Michael. Your answer looks fit better with my results.

I will ask to the sequencing company asking them about the library construction. Thank you again

Thanks you are too kind, but emergent viruses is what I do (albeit for humans not plants) and I couldn't get this one to make sense. Ask them for the Agilent library outputs, both good runs.and 50bp run. I've posted an example above.

Hi @AdriánP.L. please let us know how you get on here: one way or another we'll get this resolved if not directly via collaboration.

Dear @Michael I though It was solved, your answer looks fitting well with the trouble .

Hi @AdriánP.L. that is great news. Thanks for the update, and glad its been worked out.


text:
Best practice to perform a functional analysis on a set of VCFs from a cohort of patients
I am confused about how to perform a functional analysis after a variant calling on a cohort of patients. I have anotated all vcfs and filter them in order to get only those deleterious variants.
I was thinking on merging those vcf:
 I tried to merge them with bcftools merge and also using VariantAnnotation and transforming each sample vcf into a dataframe containing the variant lines to then row bind all those dataframes together.
As an example, Ill show some fields of teh first row of a dataframe:
CHROM   START REF ALT FILTER DP  GENEID  SIFT
  17     1000   T   A   PASS   62   TP53    0

Then, I would extract the genes affected by those variants and perform a gene set enrichment analysis, a pathway analysis similar to those performed after an RNA-seq analysis (with enrichR, for example). But I am confused about the viability of this steps. For example, I would like that, during the pathways analysis, the program I use not just reads the list of mutated genes in my cohort, but also takes into account that, if a gene is mutated in more than one sample, it should wheight more in the analysis than a gene which appears in just one sample.

As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.
This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.
Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.


text:
Outputting crunchable list of HBonds from Pymol
X-posted from the Main Stack Overflow because folks here are possibly more familiar with PyMol(?)
I have been trying to assess the strength of an interface I want to mutate using Pymol. What I am interested in is a list of contacts withing a specified interface (which I have calculated using the InterfaceResidues.py script from the PyMol wiki script library). My aim is to produce a bunch of mutants using the mutagenesis wizard and assess how the HBond network changes.
Actions>Find>Polar Contacts does the trick in visualising them, but doesn't output a file with numbers which I can then crunch to assess which mutations weaken my interface most.
I have been trying the scripts list_contacts.py and list_hb.py from Prof. Robert Campbell of Queen's University at Kingston, but with limited success (nothing gets printed or written to the file).
The scripts can be found here http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/
I am going to contact Prof. Campbell directly re: scripts but, in the meantime, and in case he's too busy to deal with a Pymol noob, does anyone have any ideas of what else I could try?
I am using Pymol 2.3.1 and Python 3.7.1

Use Biopython pdb, cross check manually in Pymol a subset of the output

@Michael it's not clear to me that Biopython.PDB would be better here than PyMOL.

How do you use these script? I had only a quick look at them: running these scripts defines a new command in PyMOL. Then you need to run that command.

Covered here: https://bioinformatics.stackexchange.com/questions/12862/extract-autodocked-protein-ligand-connections-programatically/12888#12888

@marcin I need to get the list of contacts and their distances for further analysis, not just visualise them.

Thanks a lot to all who commented/replied, the author of the scripts provided a response. It was down to a syntax change in PyMol of which I wasn't aware that caused my selection to come up empty.

Can you expound upon the syntax change, ideally giving the exact line or two of code in question? That may help others running into the same problem.


text:
isSpike function in SingleCellExperiment package is deprecated?
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use isSpike function to filter out ERCC (control) and MT (mitochondrial RNA) reads, but I am getting an error message that isSpike is a deprecated function. I am trying to use altExp instead as recommended by some people. I can't seem to figure this out. Can anyone provide some code to use instead of isSpike? I am new to bioinformatics.
I tried to swap out the isSpike function with the altExp function, but I am getting stuck on how to translate the grepl function.

I think its deprecated.. can you provide an example of your dataset? Sometimes it might have different annotation in yours

I am following the workshop here [https://scrnaseq-course.cog.sanger.ac.uk/website/index.html].
The github for the course is here [https://github.com/hemberg-lab/scRNA.seq.course].
I pulled the data from here [https://github.com/jdblischak/singleCellSeq] The data files in the tutorial are annotation.txt and molecules.txt which are listed under the Useful files section.

They hit the code at 01:59:53 in this video https://www.youtube.com/watch?v=thHgPqQpkE4&feature=emb_err_woyt

isSpike was provided by SingleCellExperiment. It is now deprecated. The notes in SingleCellExperiment indicate that one should use altExps
I am struggling with a similar problem. Can't compile a package (powsimR) because a dependency uses isSpike (not used in the code directly), so I have to hunt down where it is used by one of a very long list of dependencies.

You are on the latest version of SingleCellExperiment most likely. You should not go back to using the old version. One thing is to check out the new tutorial. In any case, the purpose of altExp() is to provide another slot of ERCC reads and MT.. etc.
So in your example (super not recommended to watch video btw) :
molecules = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/molecules.txt",stringsAsFactors=FALSE,row.names=1)

annotation = read.delim("https://raw.githubusercontent.com/jdblischak/singleCellSeq/master/data/annotation.txt",stringsAsFactors=FALSE)

We can pull out the spike in and mt genes:
library(SingleCellExperiment)
library(EnsDb.Hsapiens.v86)

Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]

Have a vector of genes to keep:
keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))

Then make the single cell object:
sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)

Use altExp to slot in the stuff:
altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

If you already have everything in a SingleCellExperiment called sce, and the old code called for something like this:
# old syntax
isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))

You can do the following instead:
# new syntax
is.spike <- grepl("^ERCC", rownames(sce))
sce <- splitAltExps(sce, ifelse(is.spike, "ERCC", "gene"))

This will remove your ERCC counts from the regular counts assay and put them in an 'alternative experiment', calling sce will print altExpNames(1): ERCC.


text:
Hypergeometic test for the overlap of two set of genomic intervals (e.g. from Chip-seq data)
Suppose I have done two Chip-seq experiments. Now after the bias correction step such as MACS2, I get two set of genomic intervals (or peaks, as what is usually called) that corresponds to each of the experiment. I want to test the significance of the overlap of these two set genomic intervals that corresponds to each experiment.
The common test employed seemed to me is the Hypergeometic test. There's even a package specifically aim for this task: ChiPpeakAnno.
Recall from the Hypergeometric distrbition -- the analogy is that we have a bag of marbles with two different colors, say one blue and another red. A sample from such distribution consists of marbles drawn from this bag, and a probability measure for this distribution is the number of blue marbles in the sample.
Could someone explain what is the correspondence between the hypergeometric test and the significance overlap between two different set of genomic intervals?

Interesting dicussion, its about sampling error in minority populations and whether this sampling error can throw a spurious "they are different result", when in fact its a sampling bias.

There is no such thing as a hypergeometric test, at least in statistical textbooks. It's a fisher test based on hypergeometric distribution.
If it is chip-seq for the same target, i.e biological replicates, significance of overlap is not quite meaningful. You get more information by for example checking the correlation of the coverage between your replicates using deeptools or anything else.
If you are comparing between two different targets, for example, chip for two different histone marks, you can use the fisher test or chi-square this way.

Partition the genome in X bp bins, these are all your marbles
Identity those that overlap with peaks from target 1 and those that don't
Do 2. for target2
construct a contingency table using 2 and 3.

So you should end up with something like this (purely made up):
                      no_overlap_target2 overlap_target2
  no_overlap_target1               9799             112
  overlap_target1                    88               1

And if there is overlap, the bottom diagonal with be higher.. less on the off-diagonal, and you can test it with a fisher.test.

So in your table, top left is the sum of the number of intervals in target1 and target2 that are not overlapping each other? Also the off-diagonal is the the number of target 1 that is not overlapping with target2 and vice versa?

sorry.. not very clear, this table is for example, we partition the genome into 10000 bins.. 9799 don't have target 1 or target 2, off diagonal is found in target2 but not in target1

Ah, I see. Thanks! Now I think I know how to do.

Yes, I have two targets. I think I'll just go ahead and use Granges package in R to do this task.

Do you mind if you could explain why this is corresponds to the hypergeometric distribution in this setting? In particular, why does "drawing the blue ball" from the urn corresponds to X-bp bins that are not in target1 nor target2?

https://en.wikipedia.org/wiki/Hypergeometric_distribution, look under working example, you are looking for the association between drawn and color

Just to make sure, I took a quick look at chippeakanno, and it uses a hypergeometric test for gene enrichment, not necessarily overlap of two chip-seq results.
As you mentioned, you can imagine in your set of genes, you have those that overlap with your peaks, and those that do not. Hypergeometric distribution calculates the probability of getting a certain number of genes that overlap with your peaks (which is the same as sampling some balls and calculate the p-value of some balls with one color). The hypergeometric test calculates a p-value to see whether this number of overlapping genes is more than expected.
I also read on Wikipedia that the P-value of hypergeometric test times 2 is exactly the P-value for Fisher's exact test. Fisher's exact test tests whether your overlap is more than (or less than) expectation. Hypergeometric tests just feel like a one-tailed version of Fisher's exact tests.


text:
Is the number of contigs and scaffolds an artifact of genome assembly and rather arbitrary, or is it canonical for each species?
E.g. bowhead whale genome has X contigs and Y scaffolds, but couldn't these easily be different numbers?

Please edit your question to describe a snippet of the whale genomes in question with examples and consider the investigation more carefully. In context the answer is it is of course arbitary. How can it possibly have taxonomic value?

This is purely arbitrary and an artifact of how well the assembly turned out. Better quality assemblies with have fewer and longer contigs.

With the ideal for your example being one contig and one scaffold or each chromosome or segment of DNA.


text:
miRNA dataset for prediction experiments
I am from the field of Computer Science and I want to perform some experiments related to miRNA and mRNA prediction.I have been searching for datasets in this area and found one called mirWalk available at: http://mirwalk.umm.uni-heidelberg.de/
when I enter, for example, the Entrez ID 10608 it appears to me the following table (I have only put the first record):

I have two questions:
a) Can I use this dataset to try to make predictions about the interaction of this gene with, for example, miRNA hsa-let-7c-5p? For what I see this interaction appears in the table, but it seems it has not been validated.
b) What are the meaning of the different fields, there is no information about it in the website of mirWalk?
Thanks for your help and consider that I am a newbie in this bioinformatics field.


text:
Get regions' information from DNA sequence data (bsgenome.hsapiens.ucsc.hg19)
I have a problem in R. I have the following dataSet (the first three rows shown) (the 5th number is the methylation level in its region):
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +

and then I'm trying to get the sequences from every row in the dataSet (427.000 rows in total):
library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet[,1]))
start<-as.matrix(as.integer(as.character(mydataSet[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet[,3]))+50)
Sequences_for_50_plus_bases<-getSeq(genome,chr,start=start,end=end)

I'd like to find a new fuction where I can have this dataSet, and find for its rows (which represent sequences in the human genome) if they refer to a promoter, or an enhancer or any other DNA regions' information, because I'd like to compare the differences between the regions' methylation levels.
Is there a way to do this?

You can get specific sequences efficiently using the following packages:
library(GenomicFeatures)
library(Biostrings)
library(tidyverse)
library(data.table)
library(dtplyr)
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

Methylation_data <- tribble(
  ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
  "chr1",15864,15866,FALSE,894,"+",
  "chr1",534241,534243,FALSE,921,"-",
  "chr1",710096,710098,FALSE,729,"+") # Arrange your dataset

Methylation_GRanges <- GRanges(  # Construct a GRanges objects
  seqnames = Methylation_data %>% pull(seqnames), # First pull in the chromosome locations
  IRanges(start = (Methylation_data %>% pull(start))-50, # then the start position minus 50 nt
          end = (Methylation_data %>% pull(end))+50), # the end position plus 50 nt
  strand = Methylation_data %>% pull(strand) # the strand your sequence is located on
)
> getSeq(Hsapiens,Methylation_GRanges)
A DNAStringSet instance of length 3
width seq
[1]   103 GCTGCTGCTTCTCCAGCTTTCGCTCCTTCATGCTGCGCAGC...CTTGGCGGATGGACTCTAGCAGAGTGGCCAGCCACCGGAG
[2]   103 ACAGTTCCAATGTAATCAGAGAGAACATCACACACACACCA...ACCTGAGCAGCACTCTGCAAAGCTGTCAAGGCGGTGAAGC
[3]   103 GATAAAATAAAGCTTAGATTGGAAAAAATATTTAAGATTCT...GGAAGCTGAGTAATTGTATGTTCAAATACTTGCAAAACAT

To write the sequences to a file, use writeXStringSet()
For the second part of what you need, GenomicFeatures has a number of functions that obtain the annotated sequences of the whole genome by its type:
transcriptsBy(x, by=c("gene", "exon", "cds"), ...)
exonsBy(x, by=c("tx", "gene"), ...)
cdsBy(x, by=c("tx", "gene"), ...)
intronsByTranscript(x, ...)
fiveUTRsByTranscript(x, ...)
threeUTRsByTranscript(x, ...)

And this is how you use one of them to find whether your sequences overlap a particular type of DNA:
Find_overlap <- function(Ref_Data,
                         Overlap_Message){

  Ref_Data <- as.data.table(Ref_Data) # Convert to data.table for efficiency

  Methylation_data <- tribble(
    ~seqnames,~start,~end,~Logical,~"Methylation Level",~strand,
    "chr1",15864,15866,FALSE,894,"+",
    "chr1",534241,534243,FALSE,921,"-",
    "chr1",710096,710098,FALSE,729,"+") %>%
    mutate(Overlap = "-",
           start = start-50,
           end = end+50) %>%
    as.data.table()

  Methylation_data[Ref_Data,
                   Overlap := Overlap_Message,
                   on = .(start >= start,
                          end <=  end)] # wherever the range overlaps, assign an identifier to it
}

Methylated_in_5UTR <- fiveUTRsByTranscript(TxDb.Hsapiens.UCSC.hg19.knownGene) %>%
  Find_overlap(., "5' UTR") %>%
  print

seqnames  start    end Logical Methylation Level strand Overlap
1:     chr1  15814  15916   FALSE               894      +       -
2:     chr1 534191 534293   FALSE               921      -       -
3:     chr1 710046 710148   FALSE               729      +  5' UTR
```


text:
What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?
What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?
Barycenter scores are calculated as 1 / (total distance from vertex v to all other vertices) in a strongly connected network. source
Closeness Centrality (Freeman) are calculated  total distance from a node v to all the other nodes in a network. source
It looks like the same centrality .
Ps:
I want to analyze closeness for my protein–protein interaction network.

First, I'll observe that Freeman centrality as you define it is the reciprocal of the Barycenter centrality as you define. It looks like people often use the reciprocal instead of the original definition, so in that sense they do look the same.
I found the following comparison here (though it doesn't name it as Freeman):

There are 2 types of distance centrality scores, Closeness Centrality
and Barycenter Centrality. Barycenter Centrality for vertex v defined
as:
1 / (total distance from v to all other vertices)
Closeness scores are calculated using the formula 1 / (average
distance from vertex v to all other vertices) and Barycenter scores
are calculated as 1 / (total distance from vertex v to all other
vertices).

So it appears that for this package it is a difference of the reciprocal of the total vs. reciprocal of the average.
Notably, wikipedia differs here, in that it matches your source:

Freeman's closeness centrality, the total geodesic distance from a given vertex to all other vertices, is the best known example.citation

Yet another source notes the existence of a "normalized" version of Freeman's closeness centrality, which normalizes centrality of a node by the degree of the network. It also comments on whether it is the reciprocal or not:

Closeness (Freeman '79)
The graph-theoretic distance of a given node to all other nodes. The sum of the rows/columns of the geodesic distance matrix of a graph.
Simple closeness is an inverse measure of centrality: the larger the numbers, the more distant an actor is, and the less central. Should really be called "farness".
Normalized version divides the minimum "farness" possible (N-1) by "farness" to simultaneously make the range 0 to 1 and invert the measure so that larger values correspond to greater centrality (truly "closeness")

So I guess that it might depend on what tools you are using. Most of the time, it might be the same. Other times, it might be normalized in different ways with either averages or for 0-1 ranges.
Closeness centrality is probably the more standard term, so it is unlikely that anyone would criticize it. But it appears that there is sometimes some variation in what people mean by that term and whether they are using a normalized or a raw closeness centrality.
In my opinion the normalization leading to a range from 0-1 is nice because it makes the measure at least somewhat generalizable across networks. But it's up to what you need to do.


text:
Speeding up network creation in WGCNA
I am following the tutorials on the WGCNA website for creating a coexpression network. It is taking way too long. Is there a way to speed up execution? All I need is to create the co-expression graph. I will then export to Python and do my analysis there.
In particular, the step "calculating module eigengenes" is taking a lot of time. I do not really care about the modules. All I want is the graph with the genes and the weighted edges.

You can use the "manual network construction" approach. Calculate TOM, then use exportNetworkToCytoscape or exportNetworkToVisANT to export the network (or simply save the TOM matrix as a table and load it in Python).

I tried doing exportNetworkToCytoscape(TOM), but I get Error in exportNetworkToCytoscape(TOM) :    Cannot determine node names: nodeNames is NULL and adjMat has no dimnames. How can you let it just use the names that were used when creating the graph?

Supply the names in the argument nodeNames.

So I did as suggested, I ran the command exportNetworkToCytoscape(TOM, nodeNames = names(datExpr), edgeFile = "edges.csv", nodeFile = "nodes.csv") but I'm curious about the result. edges.csv seems to have been populated correctly (though I'm not 100% sure, there are 5237 edges, which sounds reasonable), but nodes.csv not so much. In particular, the nodes.csv file has only 151 node entries, but I'm not sure why, since in the tutorial there are 3600 genes and so I expected 3600 entries.

If you want the entire network, you also need to set argument threshold to zero. otherwise you will only get edges whose weight are above the (default) value of threshold which is 0.5. Just be prepared to see large edge file.


text:
Weighted Network Analysis Book by Steve Horvath, taking the power of a matrix
In Horvath's Weighted network analysis, the author mentions that M^b,  where b is a real number and M is a matrix, denotes the element-wise power. This seems odd to me. The way I was taught, M^b would mean multiplying the matrix M, b times. Should element-wise power not be M^(ob) where where o denotes the Hadamard power?

You should at provide a quote since that book is not open source.

You CAN calculate M^b for b being float number, complex number or even another matrix.
As well as you can calculate exp(M), log(M), sin(M) whatever... That is very standard for math people.
Let me give two ways how you can think of it.
1) M^b = exp(b log(M)), so you need to define what is exp and log for matrices.
That can be done using power series: exp(M) = 1 + M + M^2/2! + M^3/3! + ...
log(1+M) = M - M^2/2 + M^3/3 - ...
2) More general way to calculate any function of matrix f(M) , looks like that:
if your matrix is diagonal, then it is true that: f(diag(d_i)) = diag( f(d_i)  ) .
So the question how to reduce any matrix to diagonal one.
That can be done for almost all matrices by conjugation M = C D C^{-1}.
Then f(M) = C f(D) C^{-1}.


text:
Unable to gunzip fastq files from ENA
I am trying to process fastq files in order to build gene co-expression networks (by following this tutorial). When I download any fastq file from ENA, and try to process it with the command:
ls *.fastq.gz | xargs -n1 fastqc $1, I get:
Failed to process SRR494099.fastq.gz
java.io.EOFException
    at java.base/java.util.zip.GZIPInputStream.readUByte(GZIPInputStream.java:269)
    at java.base/java.util.zip.GZIPInputStream.readUShort(GZIPInputStream.java:259)
    at java.base/java.util.zip.GZIPInputStream.readHeader(GZIPInputStream.java:165)
    at java.base/java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:80)
    at java.base/java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:92)
    at uk.ac.babraham.FastQC.Utilities.MultiMemberGZIPInputStream.<init>(MultiMemberGZIPInputStream.java:37)
    at uk.ac.babraham.FastQC.Sequence.FastQFile.<init>(FastQFile.java:80)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:106)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:62)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.processFile(OfflineRunner.java:159)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.<init>(OfflineRunner.java:121)
    at uk.ac.babraham.FastQC.FastQCApplication.main(FastQCApplication.java:316)

Similarly, when I just download the fastq files and try to decompress them using gunzip, I get, for example: gzip: SRR494099.fastq.gz: unexpected end of file
Why is this happening?

Likely that your fastq file is not successfully download based on unexpected end of file. Try zcat SRR494099.fastq.gz | tail

@PhoenixMu I get gzip: SRR494099.fastq.gz: unexpected end of file when I do this

fastqc by the way does read gzipped files, no need for decompression.

@ATpoint yeah I was just trying to see if the problem was with fastqc or with the file itself. But none of the files that I download from ENA are working

do they provide file hashes such as e.g. checksums that you can use to ensure that nothing weird happened to the files?

@MaximilianPress no but it doesn't just happen with the files I'm interested in. Any fastq file.

I just downloaded SRR494099 with the download link from sra-explorer.info and it runs fine with fastqc. Something seems to wrong with your downloads. I suggest yu download via Aspera (https://www.biostars.org/p/325010/) since this will only create the final file if download finished successfully.

Hm, that is difficult. If you can't verify the integrity of your files it will be hard to say what the issue is- but it's almost certain that they are messed up somehow. I'd follow ATpoint's suggestion or other alternative means of downloading. Possibly also switching to a clean environment with robust networking such as a fresh node on a HPC cluster might clear something up with your downloads.


text:
Global SNP frequency database?
Is NCBI's dbSNP the most comprehensive database of SNP frequency?
I'd like to get the largest possible database of allele frequencies per 'population' or geographic region. For example, on average, which country has the highest number of As at rsNNNN?
So I'm looking for global coverage above all else.

Do you mean largest in terms of number of SNPs or most individuals? gnomAD is worth checking out https://gnomad.broadinstitute.org/

Hi I agree with @sahwahn could you define the "per .... "

Looks like my answer was truncated! Sorry about that, I should have noticed.


text:
How can I convert codon coordinates to genomic position?
I am looking for a given mutation in IGV, which accepts coordinates in the form of
chr<X>:<Y>

Where X is the chromosome number and Y is the base location in that chromosome.
Eg. chr6:26091179 which has the following attributes:
Gene: HFE|LOC108783645
Marker: rs1799945
Position: chr6:26091179

Additionally, an interval can be written like
chr<X>:<Y>-<Z>

So basically from base located at Y to Z.
My questions are:
1) How do I convert HFE|LOC108783645 to chr6:26091179 (and vice-versa)?
2) How do I convert rs1799945 to chr6:26091179 (and vice-versa)?
3) On a different example, how do I convert something like codon 109 of exon 3 of beta globin gene (GTG->TTG) to chr<X>:<Y> (or chr<X>:<Y>-<Z>)?
Why are there so many ways to specify positions in genes? For me the most clear way is chr<X>:<Y>, all the others leave me guessing...

You should look at HGVS conventions used to refer to positions with respect to reference sequences. Remember, not even chr:pos notations are independent of reference sequence versions. chr6:26091179 in GRCh38 is not the same in GRCh37.

HFE|LOC108783645 is https://www.ncbi.nlm.nih.gov/gene/108783645. You can use biomaRt to go from the LOCxxxx ID to location. Searching by HFE will give you a different location though, so you might have to account for such eventualities. It is a long and often manually intensive process.
You should be able to use entrez-direct/eutils to go from rsID to position
Unless there's some sort of natural language processing API, codon 109 of exon 3 of beta globin gene cannot be parsed by existing APIs (Unless it's always codon X of exon Y of Z, in which case you could use grep too get the components. Even so, going from exon + codon + codon base change to exact base location will require a lot of wrangling.


text:
Inflated p-values in quantitative trait analysis
I am performing a quantitative trait association between the expression of one gene and ~400,000 methylation values. First, both variables are rank inverse normal transformed, adjusted for confounders (technical and for family structure) in a mixed-effects model and we use the residuals for the association. Then we do a linear regression between the residuals of expression and each methylation site, adding other four confounders in the model (age, BMI, etc.). The sample size is of ~500. For the association we are using MatrixEQTL.
When we check the p-values, they are clearly inflated. Is it something that I have to be concerned about? Did I probably miss something in the analysis? Thanks in advance.

it's a long chunk of data processing... so normally it looks like this when you are using the wrong distribution... what is the test used by MatrixEQTL ? One way to check whether that is performing ok is to permute your residuals and run MatrixEQTL again.. If it doesn't follow the diagonal, it means your data is not suitable for this.. amazing package

It performs a linear regression (least squares model) and a t-test for the significance of the coefficients.

try permuting your phenotype and see whether this qqplot still holds?

Thanks @StupidWolf

@MaximilianPress explained it better than me lol... so if you permute does this inflation go away? Honestly, I have always done it from scratch.. this chunk you mentioned above, i supposed involves some god zillion packages? Hard to tell where it goes wrong

It cannot be that your whole genome is associated with that 1 gene expression??

Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.
If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.
Some random thoughts:

What does the distribution of your methylation data look like after normalization?
Your observed distribution is basically never on the expected line, which argues strongly for confounding.
How are you adjusting for non-independence of methylation types? Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model, or at least fit 10-20 principal components to account for population structure.
How are you representing your methylation status? If it's binary you might look at this paper.

Thanks @MaximilianPress. Actually I'm don't think we have other confounding variables, we are correcting for all the possible ones, including family structure, and all data come from the same population.
After the adjustment for residuals of methylation, data is normally distributed. It is not binary, but the probability of being methylated. I'll check the papers you mention.

@SergioVillicana  You honestly can never know whether you have corrected for all confounding variables. Family structure is ok but underlying genetic relatedness, e.g. cryptic relatedness, is much more difficult to detect. Additionally, I realized that the Arabidopsis paper (mentioning the epigenetic similarity matrix) might not be appropriate, I got confused thinking about methylation as response vs. explanatory. I'm editing the answer to remove that paper.


text:
Receptor-Ligand R-package for visualization
I was wondering if there was any R-package for receptor-ligand interaction similar to SingleCellSignal (https://github.com/SCA-IRCM/SingleCellSignalR_v1/blob/master/SingleCellSignalR/vignettes/UsersGuide.Rmd)
+++ if there are vignettes and/or the graphs can be customized.
+++ if it can be linked to Seurat
Note: I have tried iTALK, but I keep getting an error that seems to be a common problem and have left a query.
Thanks

Two nice tools for receptor-ligand interaction analysis are:

CellPhoneDB by the Teichmann lab, although this is not an R package...
NicheNet by the Saeys group

Note that these tools are mainly meant for analysis of receptor-ligand interactions, although I think NicheNet also has some visualization options. CellPhoneDB needs to be run from a shell, but the results can easily be imported into R for visualization. You might check the papers on these 2 methods for some ideas of visualization. Combining Seurat's visualization options and ggplot2 should get you already pretty far.


text:
Determining a Gene's Unknown Function
I'm a beginner in this field (so pardon my rudimentary knowledge/explanation) trying to understand ways on how I could potentially identify genes of unknown function. I am trying to get as close as I can to the unknown gene's actual function.
What I've currently done, is that I've used an annotation tool to annotate my genome sample for its genes. From it, I've acquired a list of potential genes of unknown function, and I am capable of collecting its sequence. From here, I've also used BLASTn and BLASTp to align my gene sequences to find similar gene sequences.
From here, I'm absolutely clueless as to how I could proceed, since similarly sequenced genes are not always similarly functioned.
I would also like a list of other annotating tools, as well as other sequence alignment tools like BLAST, just so I could see the differences between tools.
Any help would be much appreciated.
Thank you.

Hi more detailsabout the genome, annotation and Blastp output please

Similar sequences may not always have the same function, but they're often similar enough to start with. Don't let the perfect be the enemy of the reasonably sufficient.


text:
phangorn fasta file read error
Im trying to read a multifasta sequence file into phangorn I get the following error
tipseqfile <- ("abhi_seq/clean_dup_rem.fasta")
tipseq <- read.phyDat(tipseqfile,format="fasta")

my file

Error in read.phyDat(tipseqfile, format = "fasta") :
        sequences have different length

Now is it necessary to have same length of each fasta sequence? or Im doing something wrong

You have to perform an alignment first. Phangorn use phylogenetic estimation and this requires a precise alignment position to track homologous nucleotide mutation. Maximum likelihood theory will assess the probability of a given tree for each homologous nucleotide site, if they are not alignment it can't do it and will violate the assumptions of phylogeny.
If the sequences are different lengths, which it means both nucleotide and indel ('-') sites, it means there hasn't been or there is an error with the alignment. I use muscle these days having previously used clustalo. Its really easy to use and easily configures on a linux server and its super quick. If you get stuck let me know.
I had a look at your alignment and can verify it hasn't been aligned.

Once you've run your alignment, you may additionnally think about remove poorly conserved regions, e.g. with Gblocks. It might refine your phylogeny.

damn I was using the raw fasta file..got my error I have to use the alignment file


text:
How to specifically select cluster on Featurescatter on Seurat?
Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger.
FYI, The filtered_feature_bc_matrix includes an aggregated population of tumor-infiltrating immune cells that are derived from 6 different syngeneic mouse tumors (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE121861). I downloaded fastq files and ran "cellranger count" and "cellranger aggr" commands to obtain combined data.
My strategy includes the use of the "feature scatter" function deposited in the library(Seurat). Finally, I ended up getting results below.

My question is that can we select a specific cluster(s) or identity on this plot? For example, I want to see correlation values in cluster 2 only.
The usage of Feature scatter is as follows.
FeatureScatter(
  object,
  feature1,
  feature2,
  cells = NULL,
  group.by = NULL,
  cols = NULL,
  pt.size = 1,
  shape.by = NULL,
  span = NULL,
  smooth = FALSE,
  combine = TRUE,
  slot = "data"
)
Do you know how to get the Pearson correlation of my two genes using Seurat? This is the best I can do.
Thank you.

You can use the cells argument to only display the data for a specific subset of cells.
For example to extract the values for a cluster 2:
# The cluster identity is stored in the meta.data
# In this example there is a clustering column called "RNA_snn_res.0.7"
# but any metadata column can be used

clust2.cells <- rownames(Seurat.data@meta.data[Seurat.data$RNA_snn_res.0.7 == 2, ])

FeatureScatter(Seurat.data, feature1 = "gene_A", feature2 = "gene_B", cells = clust2.cells)

Thank you so much.

You can also subset your Seurat object within the call to FeatureScatter, such as FeatureScatter(subset(Seurat.data, idents = "2"), "gene_A", "gene_B")


text:
How to incorportate RIN values as covariate in the design matrix?
I have been following the last DESeq2 pipeline to perform an RNAseq analysis with a dataset with low rin samples in the experimental (or treated) and high rin on the control ones.
I read a paper in which they perform RNAseq analysis with time-course RNA degradation and conclude that including RIN value as a covariate can mitigate some of the effects of low RIN in samples.
My question is how I should construct the design in the DESeq2 object:
~conditions+rin
~conditions*rin
~conditions:rin

none of them... :)
I cannot find proper resources where explain how to construct these models (I am new to the field...) and I recognise I crashed against a wall with these kinds of things. I would appreciate also some links to good resources to be able to understand which one is correct and why.
Thank you in advance!

Hi please edit your question to summarise the problem more clearly and include a summary of the paper in question.

hope now is more clear... :)

In the DESeq2 manual there is a section titled "How can I include a continuous covariate in the design formula?" that deals with your question.
Basically the process is no different from using a discrete covariate.
DESeq <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = ~ RIN + condition)

In this example condition is the experimental groups and RIN is the continuous covariate for the RNA-integrity.
For this to work properly there has to be a linear relationship between the RNA integrity and the gene expression values. As stated in the manual: "Continuous covariates might make sense in certain experiments, where a constant fold change might be expected for each unit of the covariate."
Therefore, it may be better to code the RIN variable as a factor like RIN_high or RIN_low. This is more robust than using the absolute values of the variable itself.
However, from your desription I get the impession that the RIN values are generally lower in the treatment group compared to the control group. If this is the case, using RIN as a covariate will by of limited use. DESeq2 will not be able to discriminate if a expression difference is due to a difference in RNA integrity or due to your treatment. Also, in this case coding RIN into high and low would result in an error because the RIN and treatment variable would define the same groups.
Generally, covariates work best when they are evenly distributed across the experimental groups.

Thank you so much for your explanation. It makes a lot of sense. I think the dataset I have to work with is going to be very difficult to safe (or at least it will need a lot of biological validation lol). I was trying to use the theory behind this paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4071332/                                                                              But they use a controlled RNA damage (with a time course) to model how it affects. I only can include rin as a covariate. I think I will go with your advice and do ~rin+condition and that's it...


text:
How do protein-protein interaction databases work?
I am a student and just trying to gain a base level understanding of PPI databases and how to interpret their data.
For example, when inputing multiple genes into String, does it only give protein-protein interactions between those identified in my gene list? Or does it compare the input gene list to a background of your selected species? If it's the former, is it normal to only want protein-protein interactions between those in your input list?
Is there a resource/text where I can find more information about PPI databases? Any help in the right direction or where I can learn more would be appreciated.


text:
Calculating polarity and hydrophobicity of a ligand/small molecule
I have lot of SMILES IDs of the ligands. What I am interested in is the polarity and hydrophobicity score of each of the ligands.
In other words, does ligand have more polar groups and can it make more polar interactions or is it more hydrophobic?
Does anybody know the way/software to look at it?
Best

You can calculate their polar surface area or logo using rdkit or openbabel. Other molecular properties may be available as well.

I see, you mean like TPSA of the molecule? Also, is there some cutoff that can be used there, in a way like TPSA of 100-150+ is considered that molecule is polar, while TPSA<100 is considered as hydrophobic molecule?

I am not sure what are the threshold or if there is a threshold or you look at literature


text:
How can I extract all known mutations of my BAM (or SNP/INDL files)?
I am using a Genome Explorer tool to see all the mutations on my own DNA.
My particular interest is on listing the variants and get their names/ids.
Here are a few screenshots:

You can see that they tell you the genomic position, the allele frequency of the mutation, the name of the disease, the Category, etc etc...
However, the filters on this tool work very badly so I cannot filter all my mutations only...
So I was trying to find out an alternative tool to get this data in such a nice way this tool provides. I tried IGV but they just provide raw readings, I didnt find a way to get a list of variants the same way this Genome Explorer does.
The closest thing I've found was Geneious:

but it crashes on my computer and is too complex to use (I wasn't able to reach that screen above, from my BAM file). And it also eats a lot of RAM that I don't have.
Are there any alternatives? Even paying?

If you have a BAM file, you could easily run a variant caller analysis software such as freebayes, though you would need a good amount of computing power like a cluster or a dedicated node with good memory. Freebayes outputs a vcf file that has all the variants listed. Then you can load it into the UCSC or Ensembl human genome browsers and look at what genes they are.

Hi! Thank you for sharing! I already have these files as well actually (.snp.vcf.gz, .indel.vcf.gz, .cnv.vcf.gz and their .tbi's). Can you detail me a bit more on the part where I can load these into some tools to list all mutations known to mankind at current 2020 (and their related disease/classification)? I am exploring both websites you just mentioned but feeling a bit lost. Thanks!

I think this one does what I want (http://www.ensembl.org/Homo_sapiens/Tools/VEP), but is giving me "bad gateway".
Any chance this is an open source project?

Found it: https://github.com/Ensembl/ensembl-vep

Yes vcf, a discussion on vcf algorithms would be cool


text:
How to remove zero value of gene on FeatureScatter plot using Seurat?
Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger. I am interested in tumor-infiltrating Tregs. I downloaded fastq files from a published paper that describes tumor-infiltrating Tregs with IL10+IL35+ (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126184), and ran "cellranger count".
Yesterday, I learned how to get a Pearson correlation value from a kind guy. I made a graph below.

In general, Tregs are defined by CD4 and Foxp3 expression. As a trial, I want to find the correlation between two genes. The correlation value was 0.05, which might have been due to the presence of Tregs that have zero value of either CD4 or Foxp3. So, my question is that can we remove cells with 0 value of genes to get more high correlation value? My data looks like this.

> Treg.data[c("Cd4", "Foxp3"), 1:30]

------------------------------------- # results

2 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGAAACCAT-1’, ‘AAACCTGAGAAACCGC-1’, ‘AAACCTGAGAAACCTA-1’ ... ]]

Cd4   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Foxp3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

------------------------------------- # results

Each dot represents 0 value.
Meanwhile, I performed to standard pre-processing workflow, as suggested in Seural Tutorial. For example,
> Treg[["percent.mt"]] <- PercentageFeatureSet(Treg, pattern = "^MT-")

> Treg <- subset(Treg, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

Can you help me? Thank you very much.

It will not improve your correlation most likely, unless the zeros are not missing at random:
library(Seurat)
data <- pbmc_small
data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^MT-")
data <- NormalizeData(data)
genes <-c("HLA-DRA","LYZ")
countdata = pbmc_small[["RNA"]]@data
countdata <- t(as.matrix(countdata)[genes,])
plot(countdata)

Now you need to subset the matrix, and remove the rows where at least one of the entries is 0, you can use rowSums(..==0)==0, meaning no zeros in the rows:
plot(countdata[rowSums(countdata==0)==0,])

cor(countdata[rowSums(countdata==0)==0,])
          HLA-DRA       LYZ
HLA-DRA 1.0000000 0.2338371
LYZ     0.2338371 1.0000000

Thank you very much, SmartWolf.
I tried to put your code on R.
I put a series of commands as follows.
(1) Tumors.data <- Read10X(data.dir = “…/filtered_feature_bc_matrix")
(2) Tumors <- CreateSeuratObject(counts = Tumors.data, project = "Tumors", min.cells = 3, min.features = 200)
(3) Tumors[["percent.mt"]] <- PercentageFeatureSet(Tumors, pattern = "^MT-")
(4) Tumors <- subset(Tumors, subset= nFeature_RNA >200 & nFeature_RNA <2500)
(5) Tumors <- NormalizeData(Tumors)
(6) genes <- c(“Cd4”, “Foxp3”)

When I put the command “countdata <- t(as.matrix(countdata)[genes,]),
it showed me that “Error in as.matrix(countdata) : object 'countdata' not found”
I tried to put another command as follows.
“countdata <-t(as.matrix(Tumors.data)[genes,])”, which at least worked with no error occurred.

sorry @raiora, was missing a line.. in your case, do countdata = Tumors[["RNA"]]@data

You are welcome. So, I am now working on it. However, in the command "countdata = pbmc_small[["RNA"]]@data", I am little bit confused. In my case, should the "countdata = Tumors[["RNA"]]@data? Sorry.

In my understanding, data refers to Tumors (seurat object) and pbmc_small would be Tumors.data

I put the command "countdata = Tumors[["RNA"]]@data", which worked.

Thank you so much @SmartWolf (Your ID does not work here for me:))

you're welcome :) sorry for the hiccup. have fun with analysis

One more question, how can I select your answer in my question? I am familiar with StackExchange site. Thank you.


text:
Selecting rows by partial match
I have a data frame one column is sample name like
Sample

TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
TCGA-2H-A9GO-01A-11D-A37B-01
TCGA-Z6-AAPN-01A-11D-A402-01
TCGA-L5-A4OP-01A-11D-A25X-01
TCGA-L5-A8NV-01A-11D-A37B-01
TCGA-R6-A8W5-01B-11D-A37B-01
TCGA-S8-A6BW-01A-11D-A31T-01
TCGA-R6-A6Y0-01B-11D-A33D-01
TCGA-R6-A8WC-01A-11D-A37B-01
TCGA-S8-A6BV-01A-21D-A31T-01
TCGA-L5-A4OJ-01A-11D-A25X-01
TCGA-LN-A7HV-01A-21D-A350-01

I want to subset samples which have partial match with these name
TCGA-IC-A6RE
TCGA-IG-A4QS
TCGA-JY-A6F8
TCGA-JY-A6FB
TCGA-L5-A43E
TCGA-L5-A4OG
TCGA-L5-A4OH
TCGA-L5-A4OJ

I have used
data[grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB", "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ", "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX", "TCGA-L5-A88T", "TCGA-Q9-A6FW"),data$Sample),]

But says that
Warning message:
In grepl(c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",  :
  argument 'pattern' has length > 1 and only the first element will be used

Likely only the first entry from the vector being grepped
Edit: As @StupidWolf suggested, I have tried:
> data[substr(data$Sample, 1, 12) %in%
    c("TCGA-IC-A6RE", "TCGA-IG-A4QS", "TCGA-JY-A6F8", "TCGA-JY-A6FB",
      "TCGA-L5-A43E", "TCGA-L5-A4OG", "TCGA-L5-A4OH", "TCGA-L5-A4OJ",
      "TCGA-L5-A4ON", "TCGA-L5-A4OS", "TCGA-L5-A4OW", "TCGA-L5-A4OX",
      "TCGA-L5-A88T", "TCGA-Q9-A6FW"),]

This works, thanks.

Answer from @stupidwolf, converted from comment:
You cannot use grepl with a vector, it must be a pattern.
I can suggest two things and you just put them together i hope, to know if something is found in a vector, you can use data$Sample %in% c("TCGA-IC-A6RE", "TCGA-IG-A4QS") etc ..
Of course data$Sample will not be in there, you need a bit of manipulation, so since your samples follow a certain pattern, you can use substr(data$Sample,....) . If you do substr(data$Sample,....) %in% c(.....) , this will give you the selection.

Not strictly bioinformatics, but a common problem:
grep assumes a regular expression, not a vector of patterns.
Assuming your patterns are stored as a vector, so e.g.
patt <- c("pattern1", "pattern2", "pattern3")

you first have to write them as a regex:
patt.regex <- paste(patt, collapse = "|")

Output of patt.regex is then "pattern1|pattern2|pattern3" which is now ready for grep or grepl:
grep(patt.regex, data$Sample)

For exact matches it would be paste(paste0("^", patt, "$"), collapse = "|")
to get "^pattern1$|^pattern2$|^pattern3$".

When the string follows certain format there is no need for partial match, see this example:
# example data
data <- read.table(text = "
Sample
TCGA-VR-A8EX-01A-11D-A36I-01
TCGA-L7-A56G-01A-21D-A25X-01
TCGA-VR-A8ER-01A-11D-A36I-01
TCGA-JY-A6FB-01A-11D-A33D-01
", header = TRUE, stringsAsFactors = FALSE)

x <- c("TCGA-VR-A8EX", "TCGA-VR-A8ER")

data[ substr(data$Sample, 1, 12) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01

Or when sample IDs lengths are not the same we can split on "-" then paste it back first 3 items, then subset as usual with %in%:
data[ sapply(strsplit(data$Sample, split = "-"),
             function(i) paste(i[1:3], collapse = "-")
             ) %in% x, "Sample", drop = FALSE]
#                         Sample
# 1 TCGA-VR-A8EX-01A-11D-A36I-01
# 3 TCGA-VR-A8ER-01A-11D-A36I-01


text:
How can i modify FASTA headers in a multi fasta file using BioPython SeqIO
I have a multi fasta file similar to this (relatively new here so uncertain of best way to present this; I have gone for an output and the code i used to make it - belt and braces...):
# create a test dataset
with open('test_dat.txt', 'w') as test_dat:
    l1 = [">AAA a aa", ">BBB b bb", ">CCC c cc"]
    l2 = ["ATGATG", "GACGAC", "TTGTTG"]

    for i in range(len(l1)):
        test_dat.write(l1[i] + "\n")
        test_dat.write(l2[i] + "\n")

# which outputs something like this
>AAA a aa
ATGATG
>BBB b bb
GACGAC
>CCC c cc
TTGTTG

I want to modify the fasta headers in my fasta file by adding some text to the beginning and, ideally, stripping out spaces (as i was advised that white space in the headings could cause problems for me down the line). My plan was to do this using BioPython and SeqIO.
I started like this:
# my attempt at altering the fasta headers using SeqIO
from Bio import SeqIO

# open my test data and create an output file
with open('test_dat.txt', 'r') as inputs, open('test_dat_out.txt', 'w') as outputs:

    # create dictionary of the fasta sequences using SeqIO
    my_dict = SeqIO.to_dict(SeqIO.parse(inputs, "fasta"))

    # renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header and remove white space
    for v in my_dict.values(): # using the values() to loop through...
        v.description = "TEST~~~TEST~~~" + v.description.replace(" ", "_")

    # write result to file
    SeqIO.write(my_dict.values(), outputs, 'fasta')

which gave me a result with some duplication in the fasta header:
#result i want
>TEST~~~TEST~~~AAA_a_aa
ATGATG

#result i get
>AAA TEST~~~TEST~~~AAA_a_aa
ATGATG

If i modify both the .id and .description in exactly the same way, as suggested here, but without removing white space:
# renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header
    for v in my_dict.values(): # using the values() to loop through...
        v.id = "TEST~~~TEST~~~" + v.id
        v.description = "TEST~~~TEST~~~" + v.description

I get almost what i want, but the header still has spaces:
>TEST~~~TEST~~~AAA a aa
ATGATG

I have tried to make the .id and .description match by adding a "_" to the end of .id, or by adding .replace(" ", "_") to both .id and .description (despite there being no spaces in .id), but I just cant seem to get exactly the output i want.
I would really appreciate any pointers on this.

You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the .id is anything up until the first white space character. One day I hope the .description and .id attributes are replaced by a single attribute like .header. Anyway, here is a solution for your case:
from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

Output on your sample input:
>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG

Much appreciated @Chris_Rands.  So my problem was that as soon as i removed all white space characters from the .description, the description - in its entirety - is considered to be the .id.


text:
Where can I visualize Ramachandran plot?
I am using Rampage webserver to obtain a Ramachandran plot. However, the server does not reveal any plot just values (favored, allowed, outlier regions).
Can you recommend other similar web sites?

Hi @polonio,  alot will understand a Ramachandran plot, I don't know anyone here who knows the webserver in question. Please either contact teh developer directly or generalise your question, by posing the problem such that it could be solved using an alternative approach, e.g. via R.

I must disagree, Rampage, albeit old is the most used site for this.

Yes, Rampage seems to be not working properly atm and I can confirm that it did show plots a few years back. Mind you the quality was not great —no plotly interactive charts etc.  @polonio210 Would you mind dropping the authors an email?

Thanks guys, I have seen recent publications and actually they work reflect Ramachandran Plots from Rampage. But, I do not know when they got them.

@MatteoFerla Yes, I will send an e-mail to them cuz the ramachandran plots of that webserver are really proffesional.

I guess the post would use a bit of clarity (what is the input, desirted output and the actualy output), but I do think it's fine here.

Here is a page with several suggested methods. I agree that Rampage looks the best of the options, but it seems that it provides some suggested workflows for drawing it in R so that it is very customizable and you can probably get it to show any features you want if you're willing to customize a little.
Here is a Python tool for drawing Ramachandran plots.
There are additionally several other servers/tools.

Thanks, it certainly appears that Rampage isnt a good option right now

Thanks, I will try using Rstudio then.


text:
about the normalization of RNAseq data for calculating distance for unsupervised learning
I have been working on clustering using RNAseq data. To compute distance, what kinds of normalization is optimal?
Can we use normalization using relative log normalization? I understand this should be used for differential gene expression analysis, not for unsupervised learning in which we do not use information of disease condition.
If you can show good papers, please share.

Hi it would be useful to describe the disease in question and some idea of the data. "Normalisation" is often used in parametric statistics to convert data to a normal distribution as a pre-requisite of the analysis, or linear modelling. Transformation is a "hit and miss" business in general. Unsupervised learning ... it depends what you mean, but you wouldn't do this for PCA or tSNE. More details please.

Typically some sort of variance stabilizing transform is used before clustering. Popular options are regularized log transformation or a vst transform, which are available in DESeq2. Note that these are NOT used for performing differential expression, just things like clustering and PCA. For differential expression one would use TMM, RLE or similar.

Thank you for your comment. Yes, I would like to calculate the distance to perform clustering, not differential expression.

So if we want to clustering after calculating distance, you mean we should use regularized log transformation or a vst transform. And when we want to perform differential gene expression analysis, we should use TMM, RLE. So we need to change normalization methods in one paper depending on analytic plan, right??

Yes, though regularized log and such are transforms and not normalizations.

So for the analytic pipeline to compute the distance, we should not use both normalization and transformation ,right?  raw count (without transformation)→vst or other transformation.
https://towardsdatascience.com/pitfalls-of-data-normalization-bf05d65f1f4c

Please see and follow the DEseq2 tutorial, this is covered there.


text:
Phylogenetics software that can represent collapsed nodes as triangles
I am trying to represent a large (100-300+ sequences) phylogenetic tree in a way that gets my point across. I figured that reducing the number of branches (by collapsing clades which include highly similar sequences) would be a good idea. The software that I'm currently using (iTOL) technically has an option for that ('collapse node'), but I'm not exactly a fan of the particular graphical representation, so I'd like to have a look at other options. Thank you for your time!

Try Figtree

For Mac Figtree here is highly recommended, it is also available for Windows and Linux here. The technique is called "collapsing clades" and is the only package I know that does it. I would recommend using it, because I would assume it is averaging the branch length of a given clade and the 'base' of the triangle is proportional to the number of taxa. Thus it recognised as being very reliable, albeit the precise algorhithm it uses ... I don't know :-O
With FigTree you get two different "triangle options" in the menu, which give different shaped triangles for collapsing the number of lineages within the select clade you chose. The big one, is the algorhithm appropriately described above. It has many options for producing publication quality diagrams, albeit I then import the output into a graphics program for further editing (there are shortcomings IMO in colouring and fonts. Then again for free is better than good.
If you are outputing a tree within a script, ETE3 (Python) will do a "collapsing of clade" option, but it isn't pretty triangle stuff. ETE3 will produce a publication quality tree diagram from this and reduce to 'low resolution', which is very useful in large analyses of loads of trees.
In summary, I don't think FigTree has any competition in the context you are seeking.
Moreover, you can configure your bin, or alias file so you can automatically open a treefile just typing,
figtree treefile.tre

Also you can reduce you font size to next point 1 or 2 (okay you can do that in anything), if you are attempting optical analysis on loads of taxa. Its also really stable. I just find it really useful, like a "jazz standard" its a "phylogeny standard".


text:
Available Protein sequence alignment dataset and HMM model
It may better to move the question here.
I am new to biology and I find my algorithm may be used in the Protein sequence alignment, since it is a henced HMM model. I find that people use HMM to generate noisy copies of the consensus sequence of different lengths. There is a figure show the process:

I am familiar with the machine learning or graphical model. To build the HMM, we need to set the model parameter, which we learn from the training data in ML. My question is whether the train&testing works for this question? If so where can I find data to training this HMM model?
It seems Professor Richard Durbin may release some dataset. But how can I find available datasets? I lost in the bio-vocabulary and fail to find datasets. Also I am wondering whether or not this topic is very minor in this community?
Update:
Since I may abuse the vocabulary, it may be better to post the original word I read.

Why don't you just build an MSA (with MAFFT/MUSCLE/CLUSTAL-OMEGA) and then build a pHMM with HMMER3? Or if you want already existing HMMs, checkout PFAM


text:
setting BLAT alignment parameters to give long matches
I'm using BLAT to align DNA sequences between similar species, and having trouble setting the parameters to return only significant matches along the 50bp query lengths.
Using -tilesize=11 -minMatch=3 -minScore=30 -maxGap=2 my understanding is that it should return alignments with at least 3 non-overlapping tiles of 11 bases, so 33 bases minimum, with gaps of up to 2 bases between tiles, and bitscore at least 30. Instead the output is showing alignments as short as 4 bases and bitscore as low as 8.
I tried setting -minMatch=5 which should not find any 55bp matches in my 50bp queries, but even that gives a handful of results, with alignment length as low as 6.
Can anyone enlighten me on the disconnect here? Thank you!
EDIT: I think I see the problem: BLAT reports on a separate line for each fragment of the match. For example, one probe has a match from 1-41 and a match from 44-50, giving scores of 77 and 13. Now I will have a different set of questions!

You're getting hits of only 5bp with BLAT, not BLAST? That should be much smaller than the word size used, how is that even possible? What BLAT implementation is this? Where can we find the explanation of what these options do exactly?

I'm running version 36. BLAT details are here: https://genome.ucsc.edu/goldenpath/help/blatSpec.html#blatUsage

Then I don't get it either. It shouldn't even be possible to get hits that short, just as you say. Could you [edit] your question and include an example query sequence so we can try and reproduce it?


text:
Why does the Genbank phylogenetic tree of 4000 sars-cov-2 sequences only display 200 nodes?
I manually selected all 4000 sequences on the 20 pages of the SARS-CoV-2 sequences on genbank here.
When I click on Build phylogenetic tree, it only has 200 nodes.
Is this expected? I was expecting to see 4000 nodes.
Is it simply a displaying issue where only 200 sequences per page are revealed while browsing the list of sequences? Is there a way to get the tree for all 4000 sequences?

Hi @shadi, I would assume the NCBI tree building program is restricted to 200 taxa/nodes. Please note a "node" is different from a "taxon" (optimal taxonomic unit).

Makes sense. Do you know if there's any option to export the full phylogenetic tree? I'm assuming that they restrict to 200 just because of the display limits.

Okay, I'll provide a formal answer, Genbank phylogeny is to get a gist of the diversity. It isn't a formal phylogenetic analysis. You can down the tree as a pdf, I don't think you can download the treefile (in fact I'm pretty certain you can't). Genbank will perform a reasonable distance method and from recollection it will "collapse clades", so if you see little triangles on the tree it represents a large number of taxa represented as a triangle.
I was surprised to learn there are 4000 SARS-CoV-2 sequences.
To make a tree you need to download the data, align (muscle) and send through a formal phylogeny program. If you think Genbank is okay, then try MegaX its very point and clicky. Generally, maximum likelihood phylogeny has returned to vogue for at least 1 decade now (previously Bayesian). However, the Beast package is a Bayesian analysis which is still going strong and the main tool now used in molecular dating.
In summary, to perform a tree analysis beyond Genbank

use the download option (there are two variants - the entire genome or the cut locus based on the query sequence)
align using muscle or clustalo (muscle is really quick(
perform a tree analysis using MEGA X, phylml, raxml or Beast depending on what you objectives are

There are 30000 sequences.


text:
Finding orthologs
I am interested in finding an ortholog for a particular gene in fungi.

Hi what are you goals for the research?

Yes, please [edit] your question and give us more context. What species is your original gene in? What fungi are you looking at? Are they close to your starting species? Far? What kind of orthology are you looking for? Do you require functional equivalence? Does your gene code for a protein and does this protein have any functional domains? You almost certainly want to start with blasting your protein against a translated genome (tBLASTn) but the details will depend on what you're after exactly.


text:
Do I need to study alternative splicing and isoformswitch separately?
I have done a study on isoform switch between different tissue type. Now, many people asking me to study alternative splicing also. But, dont they will give me same results? If there is a significant splicing difference; that will definitely leads to isoform switch. Or my understanding is incorrect ?
Thanks

Imaging isoform A uses exon 1, and isoform B uses exon 2. Now if you do an isoform switch analysis, you might see that in one condition isoform A (and exon 1) is preferred over isoform B (and exon 2). You might get similar results for alternative splicing. But if there is another isoform C that uses both exon 1 and 2, then the alternative splicing analysis of these exons are going to be affected by the expression level of isoform C as well.

Thanks @PhoenixMu