text:
Installing a package on Mac
I am trying to instll sigproextractor package on my Mac
I tried
pc-133-235:~ fi1d18$ python --version
Python 2.7.10
pc-133-235:~ fi1d18$ python3 --version
Python 3.8.2
pc-133-235:~ fi1d18$ pip --version
-bash: pip: command not found
pc-133-235:~ fi1d18$ curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py

But whatever I am trying says that can't find pip
Any help?

If you which which pip what  happens?

maybe you could try python3 -m pip install package

Sorry I am trying to install wget but says checking for GLIB... no
configure: error: Either a previously installed pkg-config or "glib-2.0 >= 2.16" could not be found. Please set GLIB_CFLAGS and GLIB_LIBS to the correct values or pass --with-internal-glib to configure to use the bundled copy.
pc-133-235:pkg-config-0.29 fi1d18$ if (!requireNamespace("BiocManager", quietly = TRUE))
-bash: !requireNamespace: event not found

At a guess you don't have python3 installed and hence pip3
https://stackoverflow.com/questions/17271319/how-do-i-install-pip-on-macos-or-os-x
Personally I would install via Anaconda3 and skip a direct Python3 installation altogether. The answer above recommends using homebrew, again personally I wouldn't touch it and go straight for a full Anaconda3 installation. You might be able to install via
conda install -c XXX sigproextractor

Where XXX is the channel, but its certainly on github. Anaconda3 is very easy to install.


text:
Translating a genome sequence into possible frames
I have a sequence below from MYC gene about which I need to translate the sequence in all possible frames AND identify (for each frame) which codons are actually used in MYC
AAAGCAATAATACAATTTAAAACCTGGGTCTCTAGAGGTGTTAGGACGTG
GTGTTGGGTAGGCGCAGGCAGGGGAAAAGGGAGGCGAGGATGTGTCCGAT
TCTCCTGGAATCGTTGACTTGGAAAAACCAGGGCGAATCTCCGCACCCAG
CCCTGACTCCCCTGCCGCGGCCGCCCTCGGGTGTCCTCGCGCCCGAGATG
CGGAGGAACTGCGAGGAGCGGGGCTCTGGGCGGTTCCAGAACAGCTGCTA
CCCTTGGTGGGGTGGCTCCGGGGGAGGTATCGCAGCGGGGTCTCTGGCGC
AGTTGCATCTCCGTATTGAGTGCGAAGGGAGGTGCCCCTATTATTATTTG
ACACCCCCCTTGTATTTATGGAGGGGTGTTAAAGCCCGCGGCTGAGCTCG
CCACTCCAGCCGGCGAGAGAAAGAAGAAAAGCTGGCAAAAGGAGTGTTGG
ACGGGGGCGGTACTGGGGGTGGGGACGGGGGCGGTGGAGAGGGAAGGTTG
GGAGGGGCTGCGGTGCCGGCGGGGGTAGGAGAGCGGCTAGGGCGCGAGTG
GGAACAGCCGCAGCGGAGGGGCCCCGGCGCGGAGCGGGGTTCACGCAGCC
GCTAGCGCCCAGGCGCCTCTCGCCTTCTCCTTCAGGTGGCGCAAAACTTT
GTGCCTTGGATTTTGGCAAATTGTTTTCCTCACCGCCACCTCCCGCGGCT
TCTTAAGGGCGCCAGGGCCGATTTCGATTCCTCTGCCGCTGCGGGGCCGA
CTCCCGGGCTTTGCGCTCCGGGCTCCCGGGGGAGCGGGGGCTCGGCGGGC
ACCAAGCCGCTGGTTCACTAAGTGCGTCTCCGAGATAGCAGGGGACTGTC
CAAAGGGGGTGAAAGGGTGCTCCCTTTATTCCCCCACCAAGACCACCCAG
CCGCTTTAGGGGATAGCTCTGCAAGGGGAGAGGTTCGGGACTGTGGCGCG
CACTGCGCGCTGCGCCAGGTTTCCGCACCAAGACCCCTTTAACTCAAGAC
TGCCTCCCGCTTTGTGTGCCCCGCTCCAGCAGCCTCCCGCGACGATGCCC
CTCAACGTTAGCTTCACCAACAGGAACTATGACCTCGACTACGACTCGGT
GCAGCCGTATTTCTACTGCGACGAGGAGGAGAACTTCTACCAGCAGCAGC
AGCAGAGCGAGCTGCAGCCCCCGGCGCCCAGCGAGGATATCTGGAAGAAA
TTCGAGCTGCTGCCCACCCCGCCCCTGTCCCCTAGCCGCCGCTCCGGGCT
CTGCTCGCCCTCCTACGTTGCGGTCACACCCTTCTCCCTTCGGGGAGACA

Where should I look for all 4 pieces of information easily? Please guide.

How do you know it's MYC gene? The answer is quite dependent on it, but in general, you probably want to just align the seuqence to a genome annotation.

The 6-frame translation  results are below you will need to do the rest. It is translated 5' to 3' but appears to contain two splice sites. You will need to use blast to resolve the start/stop of the splicing.
All 6-codon translation is available the Expasy server here. You will need to blast to resolve the splicing.

5'3' Frame 1
KAIIQFKTWVSRGVRTWCWVGAGRGKGRRGCVRFSWNR-LGKTRANLRTQP-LPCRGRPRVSSRPRCGGTARSGALGGSRTAATLGGVAPGEVSQRGLWRSCISVLSAKGGAPIII-HPPCIYGGVLKPAAELATPAGERKKKSWQKECWTGAVLGVGTGAVEREGWEGLRCRRG-ESG-GASGNSRSGGAPARSGVHAAASAQAPLAFSFRWRKTLCLGFWQIVFLTATSRGFLRAPGPISIPLPLRGRLPGFALRAPGGAGARRAPSRWFTKCVSEIAGDCPKGVKGCSLYSPTKTTQPL-GIALQGERFGTVARTARCARFPHQDPFNSRLPPALCAPLQQPPATMPLNVSFTNRNYDLDYDSVQPYFYCDEEENFYQQQQQSELQPPAPSEDIWKKFELLPTPPLSPSRRSGLCSPSYVAVTPFSLRGD
5'3' Frame 2
KQ-YNLKPGSLEVLGRGVG-AQAGEKGGEDVSDSPGIVDLEKPGRISAPSPDSPAAAALGCPRARDAEELRGAGLWAVPEQLLPLVGWLRGRYRSGVSGAVASPY-VRREVPLLLFDTPLVFMEGC-SPRLSSPLQPARERRKAGKRSVGRGRYWGWGRGRWRGKVGRGCGAGGGRRAARARVGTAAAEGPRRGAGFTQPLAPRRLSPSPSGGAKLCALDFGKLFSSPPPPAAS-GRQGRFRFLCRCGADSRALRSGLPGERGLGGHQAAGSLSASPR-QGTVQRG-KGAPFIPPPRPPSRFRG-LCKGRGSGLWRALRAAPGFRTKTPLTQDCLPLCVPRSSSLPRRCPSTLASPTGTMTSTTTRCSRISTATRRRTSTSSSSRASCSPRRPARISGRNSSCCPPRPCPLAAAPGSARPPTLRSHPSPFGET
5'3' Frame 3
SNNTI-NLGL-RC-DVVLGRRRQGKREARMCPILLESLTWKNQGESPHPALTPLPRPPSGVLAPEMRRNCEERGSGRFQNSCYPWWGGSGGGIAAGSLAQLHLRIECEGRCPYYYLTPPLYLWRGVKARG-ARHSSRREKEEKLAKGVLDGGGTGGGDGGGGEGRLGGAAVPAGVGERLGREWEQPQRRGPGAERGSRSR-RPGASRLLLQVAQNFVPWILANCFPHRHLPRLLKGARADFDSSAAAGPTPGLCAPGSRGSGGSAGTKPLVH-VRLRDSRGLSKGGERVLPLFPHQDHPAALGDSSARGEVRDCGAHCALRQVSAPRPL-LKTASRFVCPAPAASRDDAPQR-LHQQEL-PRLRLGAAVFLLRRGGELLPAAAAERAAAPGAQRGYLEEIRAAAHPAPVP-PPLRALLALLRCGHTLLPSGR
3'5' Frame 1
CLPEGRRV-PQRRRASRARSGG-GTGAGWAAARISSRYPRWAPGAAARSAAAAGRSSPPRRSRNTAAPSRSRGHSSCW-S-R-GASSREAAGAGHTKREAVLS-RGLGAETWRSAQCAPQSRTSPLAELSPKAAGWSWWGNKGSTLSPPLDSPLLSRRRT--TSGLVPAEPPLPREPGAQSPGVGPAAAEESKSALAPLRSRGRWR-GKQFAKIQGTKFCAT-RRRREAPGR-RLREPRSAPGPLRCGCSHSRPSRSPTPAGTAAPPNLPSPPPPSPPPVPPPSNTPFASFSSFSRRLEWRAQPRALTPLHKYKGGVK---GHLPSHSIRRCNCARDPAAIPPPEPPHQG-QLFWNRPEPRSSQFLRISGARTPEGGRGRGVRAGCGDSPWFFQVNDSRRIGHILASLFPCLRLPNTTS-HL-RPRF-IVLLL
3'5' Frame 2
VSPKGEGCDRNVGGRAEPGAAARGQGRGGQQLEFLPDILAGRRGLQLALLLLLVEVLLLVAVEIRLHRVVVEVIVPVGEANVEGHRRGRLLERGTQSGRQS-VKGVLVRKPGAARSARHSPEPLPLQSYPLKRLGGLGGGIKGAPFHPLWTVPCYLGDALSEPAAWCPPSPRSPGSPERKARESAPQRQRNRNRPWRP-EAAGGGGEENNLPKSKAQSFAPPEGEGERRLGASGCVNPAPRRGPSAAAVPTRALAALLPPPAPQPLPTFPLHRPRPHPQYRPRPTLLLPAFLLSLAGWSGELSRGL-HPSINTRGVSNNNRGTSLRTQYGDATAPETPLRYLPRSHPTKGSSCSGTAQSPAPRSSSASRARGHPRAAAAGESGLGAEIRPGFSKSTIPGESDTSSPPFSPACAYPTPRPNTSRDPGFKLYYCF
3'5' Frame 3
SPRREKGVTAT-EGEQSPERRLGDRGGVGSSSNFFQISSLGAGGCSSLCCCCW-KFSSSSQ-KYGCTES-SRS-FLLVKLTLRGIVAGGCWSGAHKAGGSLELKGSWCGNLAQRAVRATVPNLSPCRAIP-SGWVVLVGE-REHPFTPFGQSPAISETHLVNQRLGARRAPAPPGARSAKPGSRPRSGRGIEIGPGALKKPREVAVRKTICQNPRHKVLRHLKEKARGAWALAAA-TPLRAGAPPLRLFPLAP-PLSYPRRHRSPSQPSLSTAPVPTPSTAPVQHSFCQLFFFLSPAGVASSAAGFNTPP-IQGGCQIIIGAPPFALNTEMQLRQRPRCDTSPGATPPRVAAVLEPPRAPLLAVPPHLGREDTRGRPRQGSQGWVRRFALVFPSQRFQENRTHPRLPFPLPAPTQHHVLTPLETQVLNCIIA

thank you so much for your help. May I please ask for some more information to do this independently and confidently in future? 1. I want to stick to NCBI's tools. Is NCBI ORF finder a good tool for 6-frame translation and which parameters should I use? 2. To detect the paralogs, can you please explain the reasons behind "10,000" results and "60-80%" divergence. Also, which blast tool (blastx, tblastn etc.) should I use? 3. Can I find network of interacting proteins on KEGG ?

I've re-edited your question and my answer to ensure compliance with site rules. Long multipart questions fall into the category of "vagueness". Please read how the site works at bioinformatics.stackexchange.com/tour and
https://bioinformatics.stackexchange.com/help/how-to-ask  Ultimately it is orientated towards a third user identifying the answer from a off-site search engine, thus the title needs to be specific to the question

source: Stack Exchange

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

sort: votes


text:
Cutting the sequence into several sequences with the information of a dataframe
I have a fasta file with several 120-concatenation protein sequences. I also have a data frame with 120 names of proteins in column 1 and their length in column 2. Using this dataframe I want to create a new fasta file per sequence in which I get each protein independent with the name of the protein and the bacteria.
Here my fasta file.
Bacteria_1
Seq('GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...IEE', SingleLetterAlphabet())

Bacteria_2
Seq('GRRKDSVARVYLRPGKGEIVVNKRPVEDYFGRETLKMVLRQPLELTESLDKFDI...YK-', SingleLetterAlphabet())

Here my dataframe
      Names  LEN
0     Ribosomal_S9:  121
1     Ribosomal_S8:  129
2    Ribosomal_L10:  100
3             GrpE:  166
4           DUF150:  141
..              ...  ...
115      TIGR03632:  117
116      TIGR03654:  175
117      TIGR03723:  314
118      TIGR03725:  212
119      TIGR03953:  188

[120 rows x 2 columns]

The name of the genes and their length is the same for all my bacteria so this df is applicable for all sequences in my fasta file.
I am trying to get new fasta files like this:
>Bacteria_1 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_1 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_1 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

The second fasta file
>Bacteria_2 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_2 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_2 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

We wouldn't write the code for you from scratch (thats your job), but its an easy problem because in Python any string is an "array" or "dictionary". Thus
bacteria1 = 'GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDVIEEXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
print (bacteria1[0:120])

This will define the first gene in your list, remembering the array starts at 0 whilst your gene starts at position 1 (so you need to subtract 1). You simply parse the input list so the gene lengths form a dictionary and loop through the dictionary of gene lengths ensuring the last length is added to the current one and then issuing a write command.
You don't really need biopython in this instance, unless it will parse the gene length file you need to open (I'm not aware it can).

BTW the format you supplied isn't fasta it looks like a Biopython statement, again in this problem I wouldn't bother because Python is perfectly capable of splicing a string.

Thanks for you help Michael.

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


text:
What is the per-gene score in BUSCO output?
The BUSCO _full_table.tsv output has scores per gene (an example bellow).
In the original publication, authors just mention a condition to satisfy a bitscore, but I could not find what the bitscore actually is. In the use guide they say

If found to be complete, whether single-copy or duplicated, the BUSCO
  matches have scored within the expected range of scores ...

but again, I don't understand what the score actually represents. I am even not sure if greater the better or visa reverse. Does anybody know what the score is?
# BUSCO version is: 4.0.2
# The lineage dataset is: insecta_odb10 (Creation date: 2019-11-20, number of species: 75, number of BUSCOs: 1367)
# Busco id      Status  Sequence        Gene Start      Gene End        Score   Length
47at50557       Complete        scf1:20448-42284    20448   42284   4607.1  3554    https://www.orthodb.org/v10?query=47at50557     Immunoglobulin-like domain
547at50557      Complete        scf2:9780-23590    9780    23590   2743.8  1555    https://www.orthodb.org/v10?query=547at50557    protein unc-80 homolog
755at50557      Complete        scf3:9036-28083    9036    28083   3149.9  2206    https://www.orthodb.org/v10?query=755at50557    SRCR-like domain

Interesting question, its not immediately obvious I definitely agree

I don't think any of the guys from the Zdobnov lab are on this site yet, but do raise an issue on gitlab https://gitlab.com/ezlab/busco

i recall that the bitscore cutoffs are HMMER3 thresholds that are tailored for each individual profile (i.e. model-specific), but since i never worked on BUSCO i might be wrong

Thanks folks. If you are right about HMMER, then it is "A HMMER bit score is the log of the ratio of the sequence’s probability according to the profile (the homologyhypothesis) over the null model probability (the non-homology hypothesis)."

I might as well invite Mathieu Seppey, I think it would be healthy to have him around :-).

@Chris_Rands There you go, we got now! ^^

You are correct, it is the score of the candidate protein against the BUSCO profile hmm assessed by hmmer. If it passes a threshold, it is validated as being that BUSCO gene. The threshold is defined for each profile based on the score of the species used to produce the profile. It can be found in name_odb10/scores_cutoff. And yes, the higher the better.

Welcome, Mathieu! And thanks for the answer.


text:
Comparing phylogenetic models with different datasets
I'm a linguist interested in phylogenetic tree inference using language data. I'm posting here because I'm using Bayesian phylogenetic methods in my work (probably using BEAST and/or RevBayes). For the purposes of this question, please accept the assumption that phylogenetic trees are generally applicable to languages (I'm aware this premise may be subject to dispute). My question is intended to be methodological and not necessarily domain specific. You could imagine I'm constructing biological species trees using morphological data.
I have two datasets: A set of binary traits ($D_{b}$) and a set of continuous traits ($D_{c}$).
I'd like to compare the results of phylogenetic tree inference with (a) binary data alone vs. (b) binary data and continuous data together. Model parameters are the same, except of course (b) requires some extra parameters to model the evolution of continuous traits. In other words, both the data and model of (a) is a subset of (b).
I know this probably sounds like a strange research question. The reason it's of interest to me is that the binary data represent a data type that is well-established in historical linguistics and often used to infer language trees. The continuous data come from different language structures and haven't been used to infer language trees before. I want to know if adding this new kind of data can help aid tree inference significantly or if it has a neutral effect (or even just introduces more noise).
I'm familiar with phylogenetic model comparison using Bayes factors. But in the standard use case I'm familiar with, the dataset stays exactly the same and Bayes factors are used to evaluate the fit of different evolutionary models.
My question: Is there a standard way of evaluating the effect of adding more data to a phylogenetic tree model? Can marginal likelihoods be compared still even though the data differs in each model?
If not, would it make any sense to reframe the question as one of model comparison? I'm imagining, for example, running both (a) and (b) with identical datasets ($D_b + D_c$ for both). In (a), there would be a trivial model of trait evolution that doesn't aid tree inference at all—something like a model where traits can jump from any value to any other value with equal probability(?). In (b), there would be a more standard model of continuous trait evolution, e.g. Brownian motion.

The key issue is the data is subject to reticulation, languages are seldom clean breaks. In nucleotide data this can be identified due to deviations from physical linkage and call it recombination. Linkage doesn't apply to languages.
The fundamental assumption is one tree - the same tree - describes all the data, with reticulation that may not be a good assumption. Thus English is a hybrid language, and lots of English words find their way back into French and German 'parents'. Or perhaps a better example is Gaelic mingles with English over time. There are methods used to monitor reticulation in languages.
Thus id keep it simple. The branch length is not only directly comparable it provides linguistical information that is useful to understand. This is the amount of 'evolutionaey' change per word (it that is what you binary data is describing) and the relationship to the node ('breakpoiint') is of primary interest.
The second step is to simply draw a tree of data set A and compare it to data set A+B and data set B and the probabilities describing the sister groups. Are the trees congruent, i.e. are they statistically the same?
The area you describing is more given to parametric bootstrapping, using custom resampling, but you need to be certain that the model being used is appropriate for the data.
The final point is whether a reversible or directional matrix is being used, which impacts on the specific result.

Thanks—appreciate your perspective on this. The presence of horizontal transmission is well understood in linguistics. The general applicability of a phylogenetic tree model to language data is a big and interesting topic, but it's beyond the scope of my question here. I will edit the question to make this clearer.
It's debatable whether English qualifies as a 'hybrid' language, though certainly it has a lot of borrowed vocabulary as you point out (though again, this is beyond the scope of the question at hand).

I agree there are many interesting questions of evolutionary dynamics that can be answered by looking at branch lengths. However, it's not clear to me how looking at branch lengths will help answer my specific question (comparing the fit of a phylogeny inferred with $D_b$ vs a phylogeny inferred with $D_b + D_c$). Your second step may be one way to go, simply inferring trees with $D_b$, $D_c$ and $D_b+ D_c$ then evaluating their congruence and comparing posterior values of various nodes.

Okay its callled "horizontal transmission". In viruses it definitely isn't understood, but its calculable. Branch lengths ... the branch length and the topology are linked, hence one is an indication of the other. If the rate of evolution in one data set is disproportionate to the other, firstly you need to account for the heterogeneity in your combined data set (gamma distribution), secondly do they really give the same phylogeny?? You are assuming the rate of horizontal transfer between the data sets is approximately equal. If your data sets are big, try machine learning.


text:
Error in ruuning a package in R
I want to install Signature Tools Lib R package
I typed
> devtools::install()
Error: Could not find package root, is your working directory inside a package?
>

Packages does need this list of sotware
R package dependencies:
    VariantAnnotation,
    BSgenome.Hsapiens.UCSC.hg38,
    BSgenome.Hsapiens.1000genomes.hs37d5,
    BSgenome.Mmusculus.UCSC.mm10,
    BSgenome.Cfamiliaris.UCSC.canFam3,
    SummarizedExperiment,
    BiocGenerics,
    GenomeInfoDb,
    NMF,
    foreach,
    doParallel,
    doMC,
    lpSolve,
    ggplot2,
    methods,
    cluster,
    stats,
    NNLM,
    nnls,
    GenSA,
    gmp,
    plyr,
    RCircos,
    scales,
    GenomicRanges,
    IRanges,
    BSgenome

I have installed all of these one by one but again when trying things like the below, I get error
Error in tabToSNVcatalogue(subs = tmpSNVtab, genome.v = "hg19") :
  could not find function "tabToSNVcatalogue"
> devtools::test()
Error: Could not find package root, is your working directory inside a package?
>
> library(signature.tools.lib)
Error in library(signature.tools.lib) :
  there is no package called ‘signature.tools.lib’
>

First, we would need to know which is the package you want to install. Second this is related to Bioinformatics? it doesn't seem so, could you edit the question to include this details?

If you do ?devtools::install, you will see that the install() function from the devtools package, not surprisingly, requires a package to be specified:
install(pkg = ".", reload = TRUE, quick = FALSE, build = !quick,
  args = getOption("devtools.install.args"), quiet = FALSE,
  dependencies = NA, upgrade = "ask", build_vignettes = FALSE,
  keep_source = getOption("keep.source.pkgs"), force = FALSE, ...)

So it is absolutely normal that devtools::install() does not install anything. Just specify the packages(s) with the pkg argument.
And this is how you install Signature Tools:
install.packages("remotes")
remotes::install_github("Nik-Zainal-Group/signature.tools.lib")

Thank you but your solution gives ERROR: dependency ‘NNLM’ is not available for package ‘signature.tools.lib’

removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/signature.tools.lib’
Error: Failed to install 'signature.tools.lib' from GitHub:
(converted from warning) installation of package ‘/var/folders/wk/qwgr9yqs6ts0czx__y_w8gch0000gq/T//RtmpA8PudK/filec5bcb188c2b/signature.tools.lib_0.0.0.9000.tar.gz’ had non-zero exit status


text:
Will using smaller kmers help get larger contigs? If not, then what?
I've been using MEGAHIT to assemble metagenomes, with particular focus on specific genomic areas.
Sometimes all I get is gene fragments or pathway fragments (eg. if I know that genes A, B, C, D and E should be together, I only get A, B, C in one contig, and maybe D and E in another). That is understandable when the sequencing depth is low, and I'm assume that it's the best I can get.
However. In two of my metagenomes, the genes of interest seem to be abundant (depth of 100-300+ listed for MEGAHIT contigs). And I still have the same problem. I assume that these 'breaks' between contigs are due to natural sequence variation which breaks the assembly process. (while sequencing errors are also possible, I've used bbduk to clean my data prior to assembly)
I'd like to have larger contigs because I'm curious about the gene order. Things I've tried (did not work): visualizing assembly graphs with Bandage around Blast hits, starting with a smaller k-min (21 as opposed to 27).
So I guess I'm wondering if using even smaller k-mers might make a difference. Or if someone has another suggestion, I'm happy to listen! Thank you.
Megahit setting used:
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 1500 --presets meta-sensitive -o output

Meta-large starts from minimum kmer size 27; meta-sensitive from 21.

Typically the tradeoff is smaller kmers getting you a more connected assembly graph, while longer kmers can span repeats more and perhaps connect more contigs. So I think the logical thing might be to try and increase kmer size.

@NatWH Thank you, I can see how that would work. Considering that a MEGAHIT run goes all up to k-mer size 127 (maximum), and the results that I get are as such, shall I assume that nothing can be done just by changing k-mer size? (or even in general)

Does MEGAHIT build the results of the graph from larger kmer sizes on top of those from smaller in some way, akin to SPAdes? If so, then the method should already have accounted for this. If not, then perhaps you can try and manually increase kmer size.

You can try s-aligner (free for 15 days) which usually gets quite larger contigs and quite larger NG50 for metagenomic data containing viruses (also phages). Indeed metaSPAdes also gets larger contigs, despite shorter than s-aligner.

Ref. s-aligner: a greedy algorithm for non-greedy de novo genome assembly

Do I understand right that you are the developer of this commercial product? Please, disclose your conflict of interest.


text:
Cannot blast database created with blastdb_aliastool: BLAST Database error: No alias or index file found
I want to merge/concatenate BLAST databases as suggested here.
This is my example script that explains what I'm trying to achieve better than words can:
#!/bin/bash

for file in FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn
do
    echo "## Make DB for $file"
 makeblastdb -in $file -dbtype nucl &>/dev/null
    blastn -db $file  -query query_single_nucl.fasta -outfmt 6
    echo
done

echo "## Create testdb.nal..."
blastdb_aliastool -dblist "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -dbtype nucl -title test -out testdb

echo
echo "## Test testdb.nal..."
blastn -db testdb.nal -query query_single_nucl.fasta -outfmt 6

This is the output:
## Make DB for FAM1079.ffn
test_query  FAM1079_00001   100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM3228.ffn
test_query  FAM3228_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM6161.ffn
test_query  FAM6161_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM19036.ffn
test_query  FAM19036p_000001    100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Create testdb.nal...
Created nucleotide BLAST (alias) database testdb with 12692 sequences

## Test testdb.nal...
BLAST Database error: No alias or index file found for nucleotide database [testdb.nal] in search path [/path/to/files::]

While the individual databases work, the concatenated databases don't.

The problem was that the database must be specified without suffix, i.e. blastn -db testdb as opposed to blastn -db testdb.nal.
By the way, it produces exactly the same output as
blastn -db "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -query query_single_nucl.fasta -outfmt 6

and the performance is not better.


text:
Tools for modelling and visualizing growth of cells
I am looking for a program/GUI app/package that would help me do simulations for cell growth. I have a microscopic mathematical model, and the scenario is basically the following: I start with a rod-shaped cell characterized my its center and direction angle, it grows until it reaches a certain length, it divides into two identical cells whose angles change randomly and the process continues.
Is there any program or package that I can use for such a situation? I have found several on GitHub, but they were mostly for very complicated biological processes. I use Python/Julia/Matlab and can code, but I was wondering whether there are any tools for such simpler scenarios.

If the scenario is simpler, and you can't find anything by searching, I would consider just writing your own.
The modelling can be done in any language but if compute time is an issue - for example you want to run a simulation with many cells or many simulations - then I would use a language that can support that such as Julia or NumPy in Python. In fact Julia is a nice language for this type of thing as you end up with readable code that runs fast.
For visualisation in Julia specifically you could use the animation functionality of Plots.jl, which is basic but there, or could try Makie.jl. See the Molly.jl documentation for some examples of visualising particle-type processes using Plots.jl.


text:
Visualizing cell growth
I model the following events: a rod-shaped cell with center $x$ grows symmetrically until it reaches a maximum length/age, it divides into two identical cells, and the process continues. At division, the direction angle $\varphi$ of the daughter cells changes randomly.
So, the growth rate is $l=\frac{l_{max}}{2}+\frac{a}{a_{max}}\frac{l_{max}}{2}$
Is there any way of naively visualizing this? I found different cell modellers on GitHub but they were very complicated biologically speaking. I have some basic Python knowledge, but I have never done simulations.

Are you sure it is not 1/l ?

I would recommend looking into Python's plotting utilities, for example here.
I would recommend specifically looking at the examples that plot mathematical functions.
The basic idea is that you create a function that computes your equation based on the different parameters as inputs ($a, a_{max}, l_{max}$) and then you can plug in a range of values for each parameter and plot e.g. with $a_{max}$ on the x axis, then different lines showing a range of values for $l_{max}$, or whatever makes sense for you.
e.g.
import numpy as np
import matplotlib.pyplot as plt

def fn(a, amax, lmax):
  return lmax/2 + (a/amax) * (lmax/2)

plt.figure()
amax = np.arange(0.01, 1, .01)
plt.plot(amax, fn(.005, amax, .1), "ro") ## red points
plt.plot(amax, fn(.005, amax, .2), "go") ## green point
plt.show()

I have no idea if these parameter ranges are meaningful for you, you will have to figure that out for yourself.

Thank you for your answer! I, however, messed up in the first place. The growth rate is obviously a velocity, equal to $$s=\frac{l_{max}}{2a_{max}}$$.
After doing some reading, the length of division should be fixed, and the cells are theoretically supposed to double every 30 minutes. However, as it gets crowded, the maximum age might become bigger.

Ok, I'm not sure that I understand what you are saying but you should be able to use the tools in the same way regardless of whatever the function turns out to be.


text:
NCBI SRA database sample: control vs test
I was trying to download some data from NCBI SRA (SRA059451). There are 27 samples available for SRA059451. But i am unable to understand which samples are 'control' and 'test' samples. Please help me regarding how to differentiate the samples.

Seems to be cross-posted on Biostars, which is fine but please refer both posts to each other so that answers on one site would be utilized by the other community.

You can download the information using this link:
I did it in R:
tab = read.delim("https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=SRA059451&result=read_run&fields=run_accession,experiment_title,sra_ftp,sample_alias")

head(tab)
  run_accession
1     SRR585697
2     SRR585700
3     SRR585726
4     SRR585727
5     SRR585728
6     SRR585729
                                                                experiment_title
1 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 171
2 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 238
3 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 241
4 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 242
5 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 243
6 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 245
                                      sra_ftp
1 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585697
2 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585700
3 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585726
4 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585727
5 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585728
6 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585729
                                    sample_alias
1     Symptomatic atherosclerosis gut metagenome
2 Symptomatic atherosclerosis gut metagenome 238
3 Symptomatic atherosclerosis gut metagenome 241
4 Symptomatic atherosclerosis gut metagenome 242
5 Symptomatic atherosclerosis gut metagenome 243
6 Symptomatic atherosclerosis gut metagenome 245

Normally the information will be here, but I cannot find what does metagenome 238 etc refers to. No table in the paper that can link sample to phenotype..
Most likely you need to contact the authors. This list might be useful though.


text:
Indel quality from mpileup
Could someone please confirm that I understand this correctly? I have the res and qual columns from mpileup and I would like to match them to get the qual per base. It seems that the indel initiation (+2tt) has no qual so if I remove that I can get the phred with pythons ord().
res_with_indel = 'T*,t**,***T**t***,+2tt**,TT,.**.***.,'
res_without_indel= 'T*,t**,***T**t***,**,TT,.**.***.,'
qual = ":'21(:,838><*'915%2-;=@=99'*9'A$("
len(res_with_indel), len(res_without_indel), len(qual)

(37, 33, 33)

for base, quality in zip(res_without_indel, qual):
    print(f'{base} : {ord(quality)}')

T : 58
* : 39
, : 50
t : 49
*  : 40
* : 58
, : 44
...

To get the indel qualities then I need to match up each of the *'s across their respective indels lengths?
Is there a library that does that? Biopython maybe?
Thanks!


text:
How do I differentiate outliers from in-group variations in a DESeq2 PCA plot of 9 samples distributed into 3 conditions?
I have a PCA plot from DESeq2's plotPCA(vsd, intgroup=c("conditions")) function. I have 9 samples distributed in to 3 groups of 3 biological replicates each. My reason for trying to visualise all on the same PCA plot is to see that samples from the same group (condition) cluster together and to see a separation between the groups.
The PCA plot I have is does not appear as I expect:
1-Can I just remove a sample which I deem is the outlier and perform a reclustering of the rest on the PCA?
2-What parameter can I use to declare that a point is an outlier rather than an in-group Biological variation?
3-I also wish to find out if it is allowed to plot more than two groups of sample at a time on a PCA plot.
Thank you in advance for your kind hints.

Rather than doing any automated analysis, I like to do some exploratory analysis in these cases. Plot the counts of some individual genes from which you know that they should change and some that should not change. Also see if normalization seems appropriate on the global scale, e.g. by producing pairwise MA-plots between all the samples and see if the majority of data points is well-centered at y ~ 0 and if plots with this outliers look odd. Also try to increase the number of genes for PCA. Try to get an idea if this sample is abnormous.

Since you have so few points here, it is hard to tell which are outliers by just looking at the PCA plot.
I suggest you use correlation heatmap for your gene expression and check whether there are obvious outliers.
As for using more than one groups when plotting PCA, I think you can just try passing a vector with many groups to the argument intgroup.


text:
How to identify latent variables in single-cell RNA-Seq data
I have a single-cell RNASeq sample, in which I'd like to identify latent variables (e.g. response to stress) that I think might be affecting the clustering.
The approach I was planning to use is to reduce the dimensionality of the dataset using one algorithm among:

PCA
ICA
NMF

Then identify the most important biological processes for every component (e.g. in each Principal Component), and then finally regress out the processes that might not be relevant to the analysis or that might affect the downstream clustering.
I have a few questions:

Are there issues in this approach?
Is there a collection of gene sets (similar to MSigDB) that can be used to identify noise latent variables?
Which algorithm (PCA, ICA, NMF, or others) is more suited to this kind of analysis?

..... PCA + tSNE

@Michael why tSNE? How do you use tSNE in such case?

Its brilliant.... :-)

The larger question is how do you know you are regressing out the right stuff?

And if the so called latent variable is taking up so much of your variance... how can you be sure of the response of your interest

@StupidWolf it's a good point, however, I don't want to regress out the full Principal Component. I want to analyze which biological processes are driving the clustering and run the same analysis twice (non-regressed, and regressing only the biological processes of interest, e.g. cell cycle). My interest is more methodological (PCA, ICA and NMF work in slightly different ways, for example).

Cool.. hmm very interesting actually. I try and find the suitable scRNAseq to give it a shot. Just to clarify, yes if it is cell cycle etc we can identity it, what if it's something more ambiguous? And about regressing it out, for finding clusters etc it's not a problem. If we need to do actual inference, I would include these as covariates

I think to identify latent variables, PCA is probably not going to work. NMF might be worth trying. You might want to check out a method called consensus NMF (cNMF) (https://elifesciences.org/articles/43803) published in eLife. The author described some latent variables corresponds to cell types, while others corresponds to cell-cycle or metabolic states of cells.

Thanks! I will read it now!


text:
MaCconkey Agar colour change in non-lactose fermentors
On Macconkey agar, lactose fermenting bacteria turn the agar bright red. Why does non-lactose fermenting bacteria turn the MaCconkey agar yellow?. Does it have to do with the alkalinity of waste products?

I'm voting to close this question as off-topic because it belongs on Stack Exchange Biology. There is nothing wrong with the question, just wrong site.

McConkey agar is already red. It's trick is that it has two carbon sources. The bacteria that cannot use the more energetically useful lactose use the amino acids (peptone) generating ammonia which changes the pH and the crystal violet goes yellow. The ones that can use lactose acidify the agar a bit so go a bit more red —but not much.
This being the SE for bioinformatics, it should be said that if the genome of the organism in question is sequenced and deposited in NCBI you can find if the organism possesses the enzyme that allows the hydrolysis of lactose (beta-galactosidase) by doing a protein BLAST. The lacZ gene encoded protein in E. coli is NP_414878 in NCBI. This goes in the accession number input box. The species in question goes in the organism box. If you get something with good coverage (80%) with more than 20% homology the organism is likely lac+.


text:
what is the concept of N-channel geneChip?
I can understand two-color and one-color(Channel) gene chip; but there is another type of gene-Chip called "N-channel chips" ; I do not know how can several (> 2) samples be hybridize in one chip. I cannot find some reference for the a material for this concept. See the following code that it has two channels "E“ and "Eb"; what does that channel mean? The expression matrices seem to have the same dimension. In my understanding, the two-color or two-channel microarray is to combine two samples (treated vs. control) in a chip, and the log(ratio) is used to express the differential expression.
AE1548 = getAE("E-MTAB-1548",
               path="/Users/zhang/Documents/2020/GEOsepsis/Data",
               local = TRUE, sourcedir = "/Users/zhang/Documents/2020/GEOsepsis/Data",
               type = 'full')
AE1548raw = ae2bioc(mageFiles = AE1548)#established from raw data; multi-channel data
AE1548expr1 = assayDataElement(AE1548raw,"E")
AE1548expr2 = assayDataElement(AE1548raw,"Eb")
  AE1548expr1[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1    273812.0     97733.5      205330
2        41.0        26.5          50
3        42.5        25.0          48
4        40.0        23.5          49
5        38.5        24.0          53
  AE1548expr2[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1          25        18.0          36
2          25        18.0          37
3          25        18.0          38
4          25        19.0          37
5          24        18.5          37

Are you asking experimentally why would one design something that is not Case-Control or how technically is it achieved? The nucleotides hybridise like they normally do —there is an excess of chip/bead-linked probes to the labelled mobile targets (sample), so it does not matter too much. How you can have multiple fluorescent channels is very well documented for fluorescent dyes for flow cytometry/FACS —there litterally is a rainbow-worth of probes than can be used.

sorry for not make it clear; I update the question by including an instance. Can you help me to solve the problem? @MatteoFerla

This is a mere long comment, which may help, but I do not know the details required for the answer.
Several (> 2) samples hybridized in one chip
A traditional microarray has two parts, the mobile target, which is the labelled sample, and fixed to the filter paper, glass or doped silicon is the probe, which is not labelled and is in vast excess. You can label the target fluorescently with one of a variety of fluorophores, or with different radioisotopes, which have different energies. Nobody uses radioisotope microscopy for this except for the most obscure method, say nuclear run off assays, but I mention this just to cover all bases.
Whereas you can buy whatever fluorophore you wish, the detecting machine needs to be set up for it with lasers, beam splitters and filters. Most florescent confocal microscopes use a LED lamp, a splitter and a loads of cubes of filters. In the case of flow cytometry, you can get over 16 fluorescent channels at the same time. However, deconvolution becomes necessary when one channel "bleeds" into an other —I am not sure how one deals with it with a microarray, I suspect before the step to remove regional biases (the step with LOESS regression or similar).
Affymetrix Gene chip (R)
Now, in the case of an Affymetrix Gene chip, I do not know if there is a multichannel machine, but it could be custom build. The spots on an Affymetrix chip are 3-5 µm, so regular microscopy optics still work. Although the silicon base is not transparent like a home made slide or an Agilent microarray (if these are still made). Hopefully, someone else can answer about Affymetrix machinery!
Statistics
I am not familiar with why or how you would analyse a 3+ channel experiment. One guess could be to do a dose-response experiment. But there will be some problems:

The dye flip correction would be an issue for greater than 2 samples
Channel deconvolution may be required for dyes with overlapping spectra.

thank you for your answers. Do you know what is the two channels mean in that example? I cannot understand that two expression sets


text:
How to do BLASTP for short sequences (<20 aa) effectively?
I am trying to find sequence homology between viral sequences and my protein of interest. I have the sequences of their epitopes which varies from 5 to 500 amino acids long. For shorter sequences, it is recommended to use blastp-short but I am confused until which length it is okay to run blastp?
For shorter sequences, what are the parameters needs to be added for effective results?
blastp -task blastp-short -query POI.fasta -db virus.fasta -out output.txt -outfmt 6 -evalue 20000

You should increase the -evalue 10 default option to a higher value because this value it's adjusted to the sequence database size. Check this table, maybe this can help you to decide which parameters are the best option for you.
    OPTION      TASK          DEFAUL    DESCRIPTION
    word_size   blastp        3         Word size of initial match.
    word_size   blastp-short  2         Word size of initial match.
    gapopen     blastp        1         Cost to open a gap.
    gapopen     blastp-short  9         Cost to open a gap.
    gapextend   blastp        1         Cost to extend a gap.
    gapextend   blastp-short  1         Cost to extend a gap.
    matrix      blastp        BLOSUM62  Scoring matrix name.
    matrix      blastp-short  PAM30     Scoring matrix name.
    threshold   blastp        11        Minimum score to add a word to the BLAST lookup table.
    threshold   blastp-short  16        Minimum score to add a word to the BLAST lookup table.
    window_size blastp        40        Multiple hits window size, use 0 to specify 1-hit algorithm.
    window_size blastp-short  15        Multiple hits window size, use 0 to specify 1-hit algorithm.

Thanks! Can you please explain what is the meaning of word_size and window_size in the parameters?

The word_size parameter is used in the initial search when BLAST makes a table of all possible "words" present in the query sequence. This first search is not only for identical "words" but for all "words" that would yield an alignment with a score above a threshold, with this step, BLASTp increases speed at the expense of only a small loss in sensitivity by only trying to extend a "hit" if there is a second "hit" within a certain window_size.

How does word_size effects my BLASTp result? If I want to find 5 or 6 amino acid match at least then increasing word_size to 5 or 6 will help? It means it will try to find 5 or 6 amino acids initial match?


text:
Finding difference between groups
I have mutational catalogues of 4 samples like below
> head(out)
        Responders1 Nonresponders1 Responders2 Nonresponders2
A[C>A]A       11744          13546        2897           3655
A[C>A]C        5144           7172        2295           2163
A[C>A]G         939           1257         290            279
A[C>A]T        6065           7997        2078           2088
C[C>A]A        8969          11155        2055           2582
C[C>A]C        4050           6657        1367           1173
>

 > dim(mut)
[1] 96  4
>

Note: A[C>A]A - Single nucleotide base substitution
I want to look if mutational substitutions are different between samples
I have tried
wilcox.test(mut$Responders2,mut$Nonresponders2)

t.test(mut$Responders2,mut$Nonresponders2)

var.test(mut$Responders2,mut$Nonresponders2)
fisher.test(mut$Responders2,mut$Nonresponders2,simulate.p.value=TRUE)

kruskal.test(mut$Responders2,mut$Nonresponders2)

To find a difference but I think you may know a better test to find any difference here
This is base shift (mutation, substitution across genome of each patient). I guess mutations happen randomly but I don't know if mutational caller algorithm also pick them by randomness. I expect the shift of A (Adenine) to C (Cytosine) in a group of these patients be higher than the other group. This higher rate likely should be due to some background than randomness picked by algorithm. So I think null hypothesis is that rate of base shifting is not the same in both group. I also don't expect a normalised distribution of base shifting rather a poison distribution.
Can you help please?
Thanks

Answer from @fabio-marroni, converted from comment:
You might try to use mutationalPatterns and see if the  responders cluster together and separated from nonresponders.


text:
Possibility to save output blastn table in memory using biopython
Is there a possibility, using biopython to save output table of blastn in memory, not in file on hard drive, for process it using pandas for example and then delete it from memory?
I mean, is there a possibility to redirect output of "out" argument of this function into memory? In pandas DataFrame for example.
NcbiblastnCommandline(query="fasta1.fasta", subject="fasta2.fasta", out="output_table.csv", outfmt=6)

I always saved output tables in file before, but now I need to align one set of small fasta files against another set of small fasta files and I think that the best way is avoid writing of many small tables in hard drive because after recieving one table I have to read and process this table immediately and after that I don't need this table anymore.

You can specify send to stdout using out='-' in the Biopython wrapper.
from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd

cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()

rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', 'length',
        'mismatch', 'gapopen', 'qstart', 'qend',
        'sstart', 'send', 'evalue', 'bitscore']

data_types = {'pident': float, 'length': int, 'mismatch': int,
              'gapopen': int, 'qstart': int, 'qend': int,
              'sstart': int, 'send': int, 'evalue': float,
              'bitscore': float}

df = pd.DataFrame(rows, columns=cols).astype(data_types)


text:
Joinning columns from different files
I have 48 files in which I have
r1
Organ-specific cancer signature Contribution
Esophagus_A 2948.114
Esophagus_B 815.6483
Esophagus_C 1278.381
Esophagus_D 883.5358
Esophagus_E 1446.103
Esophagus_F 880.0416
Esophagus_G 1262.907
Esophagus_H 743.2682

r2
..
r15
and
Non_res1
Organ-specific cancer signature Contribution
Esophagus_A 2660.722
Esophagus_B 1471.017
Esophagus_C 1639.904
Esophagus_D 0
Esophagus_E 2580.775
Esophagus_F 946.3675
Esophagus_G 2013.654
Esophagus_H 0

...
Non_res33
I want a unified file in column one is Organ-specific cancer and 48 more columns are the second column of each separate file
How I can do that?

Angel! Why do you never provide example output? People ask for this in every single of your plenty of threads be it here on SE, Biostars or elsewhere. This is plain disrespectful. You seem to simply want people to do your work without investing yourself. What have you tried, why doesn't it work? You asked these kinds of simple data rearrangement questions so many times before, do you ever invest time to really understand the solutions users provide? If so you could easily find a solution yourself!

I don't want to insult any one honestly. In this question I am mentioning in output I have the common first column as the first column and second columns from each file :(

If the first column is always the same:

R cbind in a loop
bash paste file*.txt

If the first columns in not same:

R merge in a loop
bash Merge multiple files by common field


text:
Understanding BuildClusterTree of Seurat
I am trying to understand how to use BuildClusterTree of Seurat to understand the relationship between clusters. Being from neither a bioinformatics or statistical background, I am not understating how to use it. I am have the following Clustertree plot. Can someone explain how to use this function to merge cell type clusters.


text:
How to simulate phenotype from real genetic data for GWAS purpose?
I'm trying to simulate binary phenotypes from the 1000 Genome Phase 3 datasets using gcta64 --simu-cc, but no success.
Everything seems to be going well, but in the end I get:
Error: can not open the file [] to read.
An error occurs, please check the options or data

And the log shows:
Accepted options:

Here's the commands I'm using:
# Convert the VCF to plink format:
$ ./plink2 --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz --make-bed --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

# Try to simulate the phenotype:
$ ./gcta64 --bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs --simu-cc 500 500 --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

Here's the whole steps with the outputs:
$ ls
gcta64  plink2

$ ./plink2 --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz --make-bed --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
PLINK v2.00a2.3 64-bit (24 Jan 2020)           www.cog-genomics.org/plink/2.0/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log.
Options in effect:
  --make-bed
  --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
  --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz

Start time: Sun Mar 22 11:56:23 2020
16384 MiB RAM detected; reserving 8192 MiB for main workspace.
Using up to 8 compute threads.
--vcf: 220000 variants scanned.
--vcf: ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pgen +
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pvar +
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.psam written.
2504 samples (0 females, 0 males, 2504 ambiguous; 2504 founders) loaded from
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.psam.
220000 variants loaded from
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pvar.
Note: No phenotype data present.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam ... done.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim ... done.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed ... done.
End time: Sun Mar 22 11:56:28 2020

$ ls
ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed    ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam    gcta64
ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim    ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log    plink2

$ ./gcta64 --bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs --simu-cc 500 500 --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
*******************************************************************
* Genome-wide Complex Trait Analysis (GCTA)
* version 1.93.0 beta Mac
* (C) 2010-2019, The University of Queensland
* Please report bugs to Jian Yang <jian.yang@uq.edu.au>
*******************************************************************
Analysis started at 11:59:16 GMT on Sun Mar 22 2020.
Hostname: Brunos-MBP

Accepted options:
--bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs
--simu-cc 500 500
--simu-hsq 0.5
--simu-k 0.1
--simu-rep 3
--out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

Reading PLINK FAM file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam].
2504 individuals to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam].
Reading PLINK BIM file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim].
220000 SNPs to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim].
Warning: Duplicated SNP ID "rs145607083" has been changed to "rs145607083_5264"
.Warning: Duplicated SNP ID "rs145607083" has been changed to "rs145607083_5265"
.Warning: Duplicated SNP ID "rs71955229" has been changed to "rs71955229_27061"
.Warning: Duplicated SNP ID "rs71955229" has been changed to "rs71955229_27062"
.Warning: Duplicated SNP ID "rs71589472" has been changed to "rs71589472_42505"
.Warning: Duplicated SNP ID "rs563156514" has been changed to "rs563156514_49111"
.Warning: Duplicated SNP ID "rs563156514" has been changed to "rs563156514_49112"
.Warning: Duplicated SNP ID "rs539504239" has been changed to "rs539504239_79196"
.Warning: Duplicated SNP ID "rs35739849" has been changed to "rs35739849_105514"
.Warning: Duplicated SNP ID "rs148795567" has been changed to "rs148795567_123134"
.Warning: Duplicated SNP ID "rs143101359" has been changed to "rs143101359_201815"
.Reading PLINK BED file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed] in SNP-major format ...
Genotype data for 2504 individuals and 220000 SNPs to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed].
Simulation parameters:
Number of simulation replicate(s) = 3 (Default = 1)
Heritability of liability = 0.5 (Default = 0.1)
Disease prevalence = 0.1 (Default = 0.1)
Number of cases = 500
Number of controls = 500

Error: can not open the file [] to read.
An error occurs, please check the options or data

$ head ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log
*******************************************************************
* Genome-wide Complex Trait Analysis (GCTA)
* version 1.93.0 beta Mac
* (C) 2010-2019, The University of Queensland
* Please report bugs to Jian Yang <jian.yang@uq.edu.au>
*******************************************************************
Analysis started at 11:59:16 GMT on Sun Mar 22 2020.
Hostname: Brunos-MBP

Accepted options:

I'm opened to use other tools if you have any recommendation. Not sure, but looks like Plink doesn't do that (you can simulate phenotype as long as you also simulate the genetic data as well...)

PS 1: I managed to simulate it before as answered here, but I've just realised it doesn't generate reliable ratios (although the logic in place, phenotypes are still too aleatory assigned - too low precision, recall and f1 scores).
PS 2: Also asked in the BioStars forum

The question has been answered in the bioStars forum: https://www.biostars.org/p/428314/
It was missing the param --simu-causal-loci.
I created the causal.snplist file by finding out which SNP's are known causals of the treat I want to simulate.
In my case, I want to simulate a Diabetes phenotype, thus I:

Got the known Diabetes causal SNP's from snpedia
Filtered out whichcausal SNP's also exist in my genomic data:

$ cat chr3_6_10.biallelic-only.bim | grep "rs4655595\|rs6718526\|rs4402960...

Created the causal.snplist with the causal SNP's as per documentation

$ head causal.snplist
rs1111875
rs12255372
rs1801282
rs237025
...

Run gcta-simu again:

$ ./gcta64 --bfile chr3_6_10.biallelic-only --simu-cc 500 500 --simu-causal-loci causal.snplist --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out chr3_6_10.biallelic-only.gcta-simu

Got the expected results:

$ head chr3_6_10.biallelic-only.gcta-simu.phen
0 HG00096 1 1 1
0 HG00100 1 1 1
0 HG00101 2 1 2
0 HG00102 1 1 2
0 HG00103 1 1 1

The last 3 columns are phenotypes. 3 phenotypes were generated due the param --simu-rep 3.

Thanks for posting the answer to your own question, I was looking for something like this


text:
SRA run data access
Trying to download run data of Run accession SRR2155174. Here's the link: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR2155174
Please go to the link and go to the "Data Access" tab. Hereunder header "Name" you can find two links:
https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR2155174/SRR2155174.1
https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-6/SRR2155174/SRR2155174.1
Do these two links download the same data or they are for different purposes?
Thanks

It seems just like two different mirrors contains the same SRA archive. If you still doubt, you can download both (its have relatively small size) and compare their md5sum.


text:
Generate VCF from different .bam files with different chromosome names
I have two resources of .bam files. One is generated by our lab (1 sample = 1 bam). One is downloaded online (again 1 sample = 1 bam). For the downloaded samples the chromosomes are labelled: chr1, chr2, chr3 etc For our lab samples, the chromosomes are labelled: 1, 2, 3 etc. I want to generate a single VCF file of variants across all samples.
I'm using bcftools:
bcftools mpileup -Ov -f ref.fasta -b samples.txt | bcftools call -mv -o bamMge.vcf

However, I get no calls and the repeated error: [E::faidx_adjust_position] The sequence "1" was not found
I have two questions:

Is my strategy correct (i.e. bcftools mpileup)? Do I need to incorporate extra steps (NB I also have matching gVCF files)
How can I either (i) alter the chromosome labeling of one of the subsets of bam files, or, (ii) use a mapping file to match chromosome labels during the mpileup run?

Thanks

The main issue here is that your bam files have different chromosome labels. (i.e. 1,2,3 vs. chr1,chr2,chr3) as you mentioned.
This suggests that the data was mapped to a slightly different version of the reference genome. As a first layer, you might want to double check that these reference chromosomes are cross compatible (e.g. both correspond to hg19 for instance). To find more information on the reference genome used, you can check the header information using "samtools view -H "
Assuming that the reference genomes are cross compatible, you can try converting it to sam, adding the missing "chr" label in a pipe, and converting it back to a bam file. However, this will also require modification of the header file which can be challenging in itself.
The easiest (and more rigorous) approach is probably just to convert the batch of bam files with the wrong chromosome labels back into fastq files. You will then want to remap them using the exact same reference genome as the other batch of bam files. After which, you should be able to perform variant calling on the combined batch of bam files.

To add to this very correct answer, here are a couple of documentation articles that explain the key points in greater detail: https://gatk.broadinstitute.org/hc/en-us/articles/360035891071-Reference-genome and
https://gatk.broadinstitute.org/hc/en-us/articles/360035531652-FASTA-Reference-genome-format

I am not familiar with getting variants from bam files, so I will just answer the second part. From the error message, it seems that it is necessary to have consistent chromosome labels across datasets. It is also a good idea to do so because this keeps your datasets consistent and makes things a lot easier. To change the chromosomes labels, you can write your own scripts (I use Python, for example). You can also find several solutions on BioStars (https://www.biostars.org/p/98582/).

Editing BAMs or VCFs to change chromosome names is really dangerous if you're not 100% sure you know what you're doing. I've seen scripts that screwed up the file contents in subtle enough ways to mess up downstream analyses. Even if you avoid that, you can end up combining data aligned to different genome reference builds, which also ends up affecting analysis results. It's much safer to realign from scratch like KT8 describes.


text:
Per Base Sequence Content in fastqc
I have a question regarding "Per Base Sequence Content" plot for "fastqc":
In the fastqc documentation, it is written: "In a random library you would expect that there would be little to no difference between the different bases of a sequence run, so the lines in this plot should run parallel with each other."
But I don't understand why different bases in a read should follow the same pattern of allele frequency ("A/T/C/G"). I mean they are different positions in the genome and it is normal that each position has different allele.
I would appreciate it if someone could help me, please.

Heads up: This question has been cross-posted on biostars

The positions you see in the x-axis of this plot are positions on your reads, not on the reference genome. Since each read comes from a random position on the genome, the frequency of A/T/C/G on each base position should reflect the base composition of the entire genome. Therefore, the percentage of each base should be approximately constant across the read.
Hope that helps!


text:
How to compute Shanon entropy
I want to analyze some viral sequences. I want to use shanon entropy as a measure to analyze the sequences. However I do  not understand the concept of shannon entropy. Would anyone explain please how it works or advise on some resources i can use to understand shannon entropy.

There are certainly protein structural biologists on the site that could do this, however I doubt they would recommend this as a key strategy to understand the entropy of spatial relationships

I think there needs to be clarification of the question. Shannon entropy is not the same thing as the Shannon-Wiener diversity index (H), commonly referred to as the Shannon index.
The reason I brought it up is because you have the tag "population genetics". Shannon diversity index is commonly used in population genetics to compare species richness (the number of species in a given area) and their relative abundance.
Shannon entropy is more to do protein structure, which isn't really population genetics, but can relate to how stable a given mutation might be. That is not really population genetics however, population genetics is about abundance of a given allele in this context.
I am pretty sure that on this site most responders will recommend alternative free energy approaches to perform this style of analysis.

Thank you so much for your response. I totally agree with you. That was a wrong tag to place for the question i asked. I have also done some reading and realized that Shannon entropy and Shannon-Wienner diversity index are two different things. So I seeked clarification from my supervisor and he actually meant Shannon index and not Shannon entropy as he advised me before. What do i want to do; I have sequence data for viruses from different outbreaks. My plan is to divide the sequences from each outbreak into ; early, mid-phase and late phase. Then use shannon index as a measure of diversity.

Okay thanks, given the error here, it might be better to simply ask a new question. The "population genetics" was correct, but "entropy" was wrong.

How is Shannon-Wiener diversity index https://en.wikipedia.org/wiki/Diversity_index#Shannon_index any different from Shannon entropy? They are mathematically exactly the same.


text:
CRISPR Screen Control - What genes do not affect phenotype when knocked out?
I have some raw read count data from a crispr screen.
The library used was addgene's

Mouse Improved Genome-wide Knockout CRISPR Library v2

I'm following an analysis pipeline that requires counts from non target gRNA in my raw read count file.
But the gRNA library does not contain non target gRNA.
I'm wondering what genes are commonly used as phenotypic controls in CRISPR screens?
In other words, what Genes, when mutated, does not affect the cell's fitness, thus can be used as a reference when comparing gRNA counts?

I would make this in a data-driven fashion. Check those genes that are canonically not expressed based on RNA-seq data in the cell type that was used here. The respective guides should have no impact and should serve as a proxy for the technical variability of the assay.


text:
What is the difference between fixed effects and random effects in the context of Linear-mixed models?
The terms in LME i.e. fixed and random create confusion? What is the genesis that can distinguish between the two ?

A dog's response to medicine may depend  on number of days it is sick and  its sex or breed.

LME assumes two components- fixed and random. Fixed effect (mean) is constant and consistent. THE random part shows a number of variables that influence the dependent effect-size individually.

In case, we have predictor variables of continuous as well as nominal(or level-type 2 or 3) categories data, it becomes necessary to use linear as well  as non-linear (probabilistic) models. Therefore, Mixed models convey use of two types of models - linear and non-linear. Nonlinear models work under probability and linear models perform under fixed effect plus random effects.  The goof up is Mixed model.

https://stats.stackexchange.com/questions/4700/what-is-the-difference-between-fixed-effect-random-effect-and-mixed-effect-mode

Maybe worth reading the link and the article above.. I have no idea where this question is going without a context, and what is its relevance to bioinformatics?

This is a good accessible description and explanation of what random effects are https://gkhajduk.github.io/2017-03-09-mixed-models/

The definitions e.g. fixed effect  and random effects are not congruent with statistics and mathematics. Any way, thanks for posting the link.

https://www.stat.ubc.ca › BookPDF
Web results
Mixed Effects Models for Complex Data - UBC Department of Statistics - University of British ...

I do not know much about statistics but I will try my best to explain.
First, random effects are defined as the factors (categories) in the population that we are not aware of (not observed), so we are randomly sampling levels of those factors when we sample the population.
Practically speaking, random effects can be found when there are hierarchical structures in the data. For example, we want to compare the relationship between students' SAT score and their college admission rates. You might just want to run a simple linear regression. But there might be some structure in the data that is hidden. For example, the students can be grouped into classes, and then schools. You can also group students by the social-economic status of their parents.
The idea is that students belonging to different categories may have different college admission rates even if they have the same SAT score. Therefore, it is important to consider some of these "random effects" when building your model.

Linear regression hypothesis is that  effect of a
 variable is comprised as a constant component plus variability. In other  words,  the independent variable (s) is deemed to reflect a fixed (constant) part - determined as population effect size - and a variable part - random effect. Further, the term -
Linear Mixed(regression)model  is used to reflect use of continuous variables and nominal variables acting as independent variables that influence  the dependent  variable.  THE CONTINUOUS VARIABLE  is  comprised  of population effect size - and a variable part - random effect . The nominal or categories based variable effect  is probabilistic and we use the term probabilistic variable for such variables. Sometimes academics mistakenly call its effect on dependent  variable as random effect.
Finally, this may be noted that  principles of ANOVA and probability  have been combined to formulate linear Mixed effects models.


text:
Warning in fastqc
I am checking the quality control of my sequences using the Fastqc tool. For some steps, like "Per base GC content", etc., I received a warning. So, I was wondering whether I Should also take care of warning messages in the Fastqc tool or I should only care about errors!
I would appreciate it if someone could help me with this, please.

Whether the warning is heeded or not it all depends on the kind of data you are working on and the type of analysis you should be doing. If you could tell us a little bit about what you are working on and what kind of analysis you want to do that could help us

Generally yes, but it depends the type of study. Remember that Fastqc was developed for bacterial genomes, mainly, so an unusually shaped distribution in that particular dataset could indicate a contaminated library. But, for example, in an amplicon 16S amplicon dataset it's expected that other modules fail (Duplicate Sequences, Overrepresented Sequences, etc.) not only GC.

Maybe, it is good to mention that I am working on the genomic data of E. Coli. So, in this case, you believe that the warnings are also important, right?

According to the documentation "A warning is raised if the sum of the deviations from the normal distribution represents more than 15% of the reads... usually indicate a problem with the library". You need to check in the downstream analysis if you have adapters, or other kind of contamination.


text:
Since every human has a different DNA (different combinations of C, G, A, T) what does it mean to have the genome done?
I'm confused about the difference between genome and DNA. Is it correct to say that the same type of bacteria has the same DNA? But my understanding is that it is not correct to say that the same type of human has the same DNA, since every human has a different DNA. What am I missing here?

I'm voting to close this question as off-topic because it is a basic biology question and is not really bioinformatics. Biology SE would be a better place to ask this.

Closing is never a good idea, unless it is something outrageously wrong or off-topic. Not the case here.

As the OP, you'd surely have a bias against closing your own question.

@RamRS, It's clear to me that the OP is asking what people mean when they say "the honey bee genome," "the horse genome," or "the human genome," given that most individual honey bees, horses, and humans differ genetically. Although I'm neither a biologist or a bioinformaticist, I have the impression that genome projects were deeply intertwined with the rise of bioinformatics as a field, so this looks like a bioinformatics question to me. It may be a basic one, but to me it feels basic in the sense of being fundamental.

Seems fairly obvious to me: they got a human genome.

@user253751 The current top rated answer reaches a different conclusion.

I work in genomics but am not a geneticist, so I also had this question for many years. The answers below are correct -- the human genome is a reference point from which we discover variation in others -- but it underlies how confusing genomics jargon can be. Another example is "variant calling". Why "calling"? I was told it was the meaning of "to call" as in "to call a race" and even then it took me a while to understand! Jargon is coined by metaphor, so in my opinion, genomicists must have quite an "interesting" view of the world!

All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are also often fatal, so the random mutations we see of living adults are much more restrictive than all the random mutations that occur.)
To identify the human genome is to study many people's DNA and to label the parts that are shared between everyone, the parts with two or three variations, and the parts with even more variations. Even though every person has different DNA, we can still say they fit the pattern, just like every T-shirt may be unique but they all fit the T-shirt pattern and not the trousers pattern.

DNA is only the substance that makes up the genome.
DNA can come with many forms. For example, the fragments you get by PCR in the lab are DNA. A bacterium also has DNA molecules in it. You may also find DNA fragments in many other places, like in our blood, or at a crime scene.
Genome, however, is a very specific term. It means all the heritable substances in an individual or an organism. Most times, when we talk about the genome in eukaryotes, we refer to the nucleus genome, that is, the DNA in the nucleus.
Note that genomes do not have to be made of DNA. RNA viruses, for example, only have RNA as their genomes. Yet we still refer to that as a genome, because that is the heritable substances for those viruses.

DNA is the product you get if you paste the nucleotides A, T, C and G together. This is DNA and all genomes consist of these basic building blocks like a book consists of letters. Still, the genome, or the content of the book is the specific combination of these letters plus the proteins that bind to the chromosomes, the folding of the DNA and all the epigenetic modifications that can happen. So yes, we all have the same DNA in terms of the same building blocks, but the content is different.

Might have misunderstood you - but "we" have sampled and "decoded" the DNA of a few humans (often volunteers who surrendered their body to science after their death) - but this gives of course only the genetic code of those few individuals.
When scientists say we decoded the genome they mean they deciphered which parts of the DNA "do" specific things - the actual "genes" or gene-sequences - usually done when looking at the DNA molecule packed up into the chromosomes to give easier "orientation".
And by my understanding the first task was accomplished - supported by use of supercomputers, but the deciphering of the gene-sequences still has a long way to go for scientists to state "complete." But there is growing success - And as usual scientists will not stop until complete.

Why do you think we have only/mostly sampled the DNA of the deceased? You can take a swab of your DNA (non-harmful, not even painful) and send it to various labs for identification of inherited diseases, genealogy (ancestry), etc.

I wrote about those as I wanted kinda honor those who donated their bodies - but you are right .. of course the DNA samples can come from living humans too since literally a tiny drop of blood would be enough @CJDennis

That makes you a good human being, but it's not good to have emotion in answers. Answers are best when they stick to the facts.

"To have the genome done" means that we have sequenced all the DNA in particular reference individuals. You are correct that every human has different DNA, so there is no such thing as a single "human genome", as (almost) everyone's is unique. But the vast majority of everyone's DNA is identical - that's what makes us human. Having the full sequence of a few "normal" individuals who don't have glaring genetic defects is extremely useful for learning about the genome of all individuals, since so much of it is the same from person to person.
I've used the term "reference" a few times, since we use a few individuals' genomes as an expectation of what a human genome looks like. The process of sequencing DNA involves breaking it into many bits, and then recombining them using a computer. You can put the pieces together without the reference, but it is much more difficult and time consuming. Even if there are differences from the reference, having an expectation of what >99% of the genome should be makes the problem much simpler. In a way, it's like the difference between doing a jigsaw puzzle without knowing what it's supposed to look like, and using the picture on the box as a reference. Even if your exact jigsaw puzzle has a few different pieces, the box reference is  mostly the same and immensely useful.

Think of it as a special book. Everyone's book is different, but they all have the same number of chapters and each chapter covers the same topic in each individual book (sorta). The problem is we can only read a few words at a time.
When the genome is "done", it means that someone has read a bunch of these sentence fragments from one person's book and then assembled all the words in the correct order to make a complete book. Then we can get a good idea of how many chapters it has and what those chapters are about.
However, the human genome wouldn't be a novel, it would be more like a set of encyclopedias.

When people say "the genome" for a species, they mean a reference genome. To quote wikipedia,

A reference genome (also known as a reference assembly) is a digital nucleic acid sequence database, assembled by scientists as a representative example of a species' set of genes. As they are often assembled from the sequencing of DNA from a number of donors, reference genomes do not accurately represent the set of genes of any single person.

The genome is updated over time as more is learned, the most recent human reference is hg38.
For individualized results, the patient is sequenced and the data is compared to the reference genome. Then, the differences are presented as a list of genetic variants. This is more efficient than looking at the entire genome for each patient, because on average greater than 99.9% of sites match the reference.

Yes you are right.  99.9% of our genome remains same, it is only the 0.1% which makes each individual unique from the rest of us. Single nucleotide polymorphisms (SNPs) contributes to the 0.1% variation in our genome. SNPs explain the variation in eye color, susceptibility towards disease conditions etc


text:
Is there a computational tool or possibility to identify mRNA isoforms from the count matrix of a bulk RNA sequencing dataset?
I have the counts matrix of an RNA sequencing dataset of fibroblasts and I wish to identify isoforms of a particular gene of interest in it. Can anyone please hint me on a bioinformatics method to achieve this?
In literature, 9 different isoforms of my gene of interest have been reported. My colleague who does the wet lab experiments has identified 3 of those isoforms through western blotting in the fibroblasts that later got sequenced. There might be more isoforms in these fibroblasts.
I will be grateful for any computational hint to identify gene isoforms from my dataset, irrespective if they are differentially expressed or not.
Thank you in advance for your kind help.

What tool did you use for quantification? Most tools produce transcript level counts, so this should be a pretty straightforward task once you get a list of the identifiers of the transcripts you're interested in.

If you have only counts on a gene level, you might need to go back to the bam file to get the isoform quantification.. Something quick might be kallisto?

You can back to your bam file and use cuffdiff tools it is very able to do your job follow this to liks
http://cole-trapnell-lab.github.io/cufflinks/manual/
https://github.com/Jeanielmj/bioinformatics-workshop/wiki/Differential-Analysis-with-Cuffdiff
My best

As a remark, since it is 2020 you better you a recent and non-deprecated tool such as salmon or kallisto to get transcript level abundance estimates and then do differential transcript expression analysis with tools such as DEXSeq, swift or sleuth (the latter is downstream of kallisto from the same group).

If literally all you have is the count info by gene, you can't magic up the isoform breakdown.


text:
Location of Reads in Sequencing
I have questions regarding sequencing:
Are paired-end reads from different stands?
What about single reads, are they coming from both strands?

Mates in a pair are giving signal for a single strand, though this is only apparent in bisulfite-treated samples. The reason for this is that you are sequencing the ends of fragments of material and it's random which strand of this (or both in different clusters) will get sequenced. Different pairs of reads will then randomly sequence one strand or the other.
Single-end reads are the same as paired-end reads, the second read in the pair is just not made.

Now, I understood from you that the mates in paired-end reads are coming from two different strands, but based on the figure in https://www.biostars.org/p/95803/ and also some other sources, I implied that they are coming from two different stands as they are read from two different directions. Could you please help me in understanding this?


text:
Can I automate CLUSTALO and output alignment sequence identity?
I've detected homology between targets of ligands in drugbank and proteins in the proteome of a pathogen. I've parsed the output very rudimentary and calculated my query coverage. This exists in an excel file (see attached)
I'd like to further filter these hits. For query coverages > 70% alone I have 1500+ hits. Some redundancy exists and some hits are for drugbank targets for which no data has been deposited in the PDB.
So there are ways of further filtering these hits, but I'd like to do pairwise alignments first before I do anything else, as I think this is the most rational filter I can think of.
I do not want to filter based on lack of PDB data necessarily, as when I check my high query coverage hits, I'm getting predictions that are actually backed up by in vitro data. This includes high coverage hits lacking x-ray data
My question:
I have CLUSTALO installed on my desktop: How do programmatically send it jobs if I have the Drugbank sequences and the pathogen sequences in the same file? I would also like to parse out the pathogen IDs from the excel file. There’s 1500 and they are a bit scattered throughout the columns.
0.980634         -                    572         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-263             874.5          0.0            1        1      6.6e-267  1.4e-263   874.2      0.0       5         562    15     570    12     572    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005678  PE=3             SV=1

0.950704         -                    580         sp|P06672|PDC_ZYMMO                         -                     568                               1.2e-258             858.0          0.0            1        1      6.8e-262  1.5e-258   857.7      0.0       6         546    10     560    7      578    0.98   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_005014  PE=3             SV=1

0.950704         -                    691         sp|P06672|PDC_ZYMMO                         -                     568                               8.2e-162             538.3          0.1            1        1      2.4e-164  5.1e-161   535.6      0.1       1         541    89     650    89     658    0.89   Acetolactate                            synthase                  OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_006830  PE=3             SV=1

0.492958         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            1        2      3.3e-83   7e-80      267.8      0.0       3         283    2      287    1      292    0.94   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

0.107394         -                    585         sp|P06672|PDC_ZYMMO                         -                     568                               8e-87                290.7          0.0            2        2      1.7e-08   3.6e-05    21.1       0.0       413       474    436    505    376    583    0.72   Uncharacterized                         protein                   OS=Lomentospora            prolificans        OX=41688             GN=jhhlp_001897  PE=3             SV=1

you will need to structure your questions so it is one question per post by breaking-down your goals into components. At present your question is not answerable I'm afraid, the site rules are designed to enable search engines to easily identify the post.

Is this an automated comment? I shall parse this out into multiple questions (pun intended)

Thanks @Mike (thats not automated "thanks")

Also, make sure to tell us what operating system you are using, and why excel (!?) is involved. Finally, explain if you need every possible combination of pairwise alignments from your sequence file or only some.

Thank you, Michael! Terdon, I just wanted to view the output in a consolidated table so I could send to others for review, though of course the awk one-liners and GREP commands were done on the command line (gnome emulator); I'm using the latest version of ubuntu

I'm using CLUSTALO mainly for single alignments; I have a few MSAs I need to run, but not many. Sorry for the confusion. I'm using CLUSTALO for the single alignments we've agreed it's gotten us optimal alignments

If you can give us a few entries from your input files and explain what you want to do with them, we can try and answer. But I don't understand why you would be using an excel file. Can you export to a simple text-based format? What do you have in the excel? We can't really help you process data if you don't show us the data you need to process.

Sure thing, see an edit to the original post. The original output is in a .txt, I just copied and pasted it into the excel sheet to send to a colleague.

Peronsally although I automate alot and automate clustalo I don't automate downloads because Blast online is pretty slick and has very neat tools, but its context specific.

Sure understood, automating clustalo, which I have locally, will do me well right now. I would eventually like to look at automating what looks to be a simple pipeline

So, as the first reply said, there's sort of multiple questions here.
For filtering things out from a table, I would absolutely use pandas. https://swcarpentry.github.io/python-novice-gapminder/08-data-frames/ Pandas is a way of managing dataframes in python. It looks like what you want from your table is to remove certain rows, if the content of those rows at a certain column matches certain IDs from a list (pathogen IDs). That's the sort of thing the pandas library is good at.
For managing workflows, I use either bash or snakemake. Snakemake is a python thing. It has some advantages, such as being able to start from an intermediary file (for example, if a workflow is "sequences -> multiple alignment -> phylogenetic tree", it can easily start from the multiple alignment). Also deals with branched workflows quite well. Can be easily installed with conda (anaconda is absolutely a lifesaver for bioinformatics work imo). I still use bash for simple workflows, because it's slightly simpler/faster.
For potential filtering - If you want to keep the top hits as they are, but filter the lower-scoring hits based on whether they have a PDB ID or not, I'd personally look into Entrez queries. They're a way of making your search more specific. Here they're shown for WebBlast, but they work quite well for command line blast too. https://www.youtube.com/watch?v=bxx5uaKjMa8

Thanks so much for this; my issue is, or my question is, will entrez have the CDS-based proteome that I'm using? I have a proteome from deposited into EMBL-EBI

@Mike I'm afraid that I haven't used EMBL-EBI much, so I'm not sure how to help. Entrez blast queries are mostly an NCBI thing. Still, the 'sequences associated with solved protein structures' should be the same in both cases, because they represent links to the PDB protein database. So I guess there are different options, from re-doing the search in NCBI, to downloading the sequences of interest from EMBL-EBI and then blasting in the NCBI with Entrez (either web-based or command line), to trying to get the NCBI code for EMBL sequences...

Thank you! I swear I replied to this, so sorry Laura

The question is still in numerous parts,

parsing a blast output
obtaining sequences
feeding this as input into clustalo

If I was automating this I would parse the Blast genbank codes [regex either Perl one-liner (scripts aren't trendy ;-( ), or Pythons "re"] and feed them into NCBI's efetch to hook out the sequence data. Then I would create a subprocess that feed the files into clustalo. In Perl you simply issue the linux command wrapped in "`" either side, in Python use
import os
os.system('clustalo -in file.in -out file.out')

Ensuing your 'bin' is correctly configured to execute clustalo
or
import subprocess # use popen

To my knowledge there is no wrapper for clustalo unlike several other phylogenetics programs in BioPython, which enormously simplifies pipelines. Others might recommend using shell rather than e.g. Python.
BTW I personally use Muscle having been a long term user of clustalo, thats a long discussion.

To be honest you should look at the power of NBCI's online blast utility. The download facilities are certainly powerful particularly if you are pulling down e.g. 1500 sequences from blast hits, because it automatically can splice the hits of interest from their genes/genomes. However, its context specific.

To be clear, I already have the sequence data for these uniprot IDs on my desktop in another file and they are linked to their sequences. I wish to do pairwise alignments on sequences detected by jackhmmer

Hmmmmm  ... anyway subprocessing is the specific answer you requested on how to automate clustalo ... unless you are using the bash shell

Currently doing some bipython tutorials at the moment, but yeah I do a lot through the command line

Biopython does include a wrapper for ClustalO (and MAFFT, MUSCLE etc.)


text:
Overlap in Paired-end Reads for Sequencing?
Regarding Sequencing:
Do/can the paired-end reads have overlaps?

Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.


text:
Confusion about GISTIC threshold
Please, can you help me with a confusion?
In whole genome sequencing, if
Total copy number (major + minor allele copy number by SCAT R package) for SMAD4 in a sample = 4
Th estimated ploidy by SCAT for this sample = 4.043696
So, the segment mean (required for GISTIC) for SMAD4 for this given position of genome in the given sample = log2(4/4.043696) =  -0.01567461
We can consider this as a loss of copy number for SMAD4
But with GISTIC threshold -td = 0.1 this is neutral copy number change (if I am not wrong)
So, If I want to ask GISTIC to give me
amplification: copy number ≥2 × average ploidy
gain:  copy number >1.25 × average ploidy
loss: copy number <0.75 ×average ploidy
deletion: copy number 0

If I put -ta = 1.25 and -td = 0.75 , do I achieve the above?
I am seriously struggling with this intuition of relationship between segment mean and GISTIC copy number threshold
Thank you for any help

You've invested quite the effort in this question, Angel. I hope this gets answered soon and answered well. Cheers!


text:
Filtering genes from cuffdiff results
This question has also been asked on Biostars
I have run cuffdiff (with statistics turned ON) to compare two groups of samples: Control group and Late AD group.
This is the command I ran to be precise:
cuffdiff -L Control,AD_Late_Braak -p 8 --total-hits-norm --frag-bias-correct ../References/ensembl.GRCh38.99.fa --multi-read-correct --library-norm-method quartile ../References/Homo_sapiens.GRCh38.99.chr.gtf Early_Braak_Control_1,Early_Braak_Control_2,Early_Braak_Control_3 Late_Braak_Sample_1,Late_Braak_Sample_2,Late_Braak_Sample_3

I'm looking at the output from cuffdiff and I see a gene_exp.diff file which contains the results of the differential expression testing. I want to know what is the best way to filter the results of this gene_exp.diff file so as to restrict the number of genes that are up-regulated and down-regulated to a list between ~50-850.
P.S My thoughts are to tweak the p-values and log2 fold change values, but seems like a trail-and-error method, so I was wondering if there was a more "formal" method/approach?

To filter your gene_exp.diff either you export it to excel and then filter because it is columns separated by tabs or use the command (awk and head -n) to filter and have the number of lines you want if you master Linux


text:
Marking optical or PCR duplicates with picard vs. samtools flagstat
I am trying to identify the best library prep method for noninvasive prenatal test samples, to be sequenced with Illumina Novaseq.
One metric that I am evaluating is the number/ percentage of PCR or optical duplicates. In order to achieve that I have done the following steps:

aligned the fastq reads against the reference genome with bwa mem
the .sam output from bwa mem was converted to bam, sorted and indexed
the duplicates were marked with the command MarkDuplicates from picard

Then if I call samtools flagstat on the sorted bam file which had the duplicates marked with picard I get:
26595942 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
466 + 0 supplementary
1636809 + 0 duplicates
24969064 + 0 mapped (93.88% : N/A)
26595476 + 0 paired in sequencing
13297738 + 0 read1
13297738 + 0 read2
21102678 + 0 properly paired (79.35% : N/A)
24049752 + 0 with itself and mate mapped
918846 + 0 singletons (3.45% : N/A)
462212 + 0 with mate mapped to a different chr
182138 + 0 with mate mapped to a different chr (mapQ>=5)

Which sounds a bit strange to me, as no duplicates were found.
However, the MarkDuplicates command also outputted a metrics file (dups/dupsMetrics.txt) that did identify some duplicates. This was achieved with the following command:
java -jar picard.jar MarkDuplicates I=alignment/sample.sorted.bam O=alignment/sample_markDup.bam M=dups/dupsMetrics.txt

I see that the value under the column PERCENT_DUPLICATION in dups/dupsMetrics.txt is 0.065555 (given the column name, I would expect that the value is indeed 0.06%, and not 6.55%)
What would be the explanation between the discrepancy in the output of samtools flagstat and MarkDuplicates. Is it reasonable to expect that the value displayed by samtools flagstat is rounded to zero, given the low percentage which was observed by MarkDuplicates (0.06%)?

There are duplicates, in this line:
1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356
According to samtools documentation for flagstat:

Provides counts for each of 13 categories based primarily on bit flags
  in the FLAG field. Each category in the output is broken down into QC
  pass and QC fail. In the default output format, these are presented as
  "#PASS + #FAIL" followed by a description of the category.

So you can usually take whatever is in the first column

@StupidWolf's answer is correct -- that first number in the flagstat output is what you want to look at to see the number of reads marked as duplicates. I wanted to add that the number given in the Picard metrics file is in fact 6%, not 0.06%. This is due to an ambiguity that is widespread in Picard metrics files; in many places the program emits a fraction even when the column name is "percent of something". Soo Hee Lee wrote a very short GATK blog post about this way back in the day, see the archived version here.

Good to know that the percentage of duplicated reads from picard and from samtools match