text:
Huge variability in blastn speed between same short sequence job
I have been running a couple of blast searches and I am very suprised and somewhat frustrated by the huge variability in job completion time. The specifics of my blast search are as follows:

Running blast using webbrowser (chrome) via ncbi
19 nucleotide query
database: WGS of one particular sequencing project/bacterial strain (so database is VERY small, in total approx 4.5 MB divided over 60 contigs)
parameters are mostly untouched and are automatically optimized for short queries by blast

In one instance, I started the job more than an hour ago and it is still not done. Knowing from experience that this previously hasn't taken that much time, I started a new job in a different browser tab with exactly the same sequence, database and parameters, and now it finished within 20 seconds. What on earth is happening here?

Welcome to Bioinformatics.SE! I assume you're running BLAST on the command line, and not using a web browser to connect to NCBI's BLAST website. Are you running this on the computer on your desk, or are you connecting to a remote computer? It's possible the machine you're connecting to is experiencing intermittent disk issues.

You mention "tabs" so what service are you using? Please [edit] your question and tell us exactly what you are running and how. What blast? If you're using a website, which one? If not, what blast executable? On what kind of computer?

Thanks for the tips, I edited my question a bit providing more information - running blast via NCBI (in both cases)

First of all, if you're using an online service like NCBI's BLAST, then the time to completion will depend on the load the service is experiencing. Your job is ut in a queue and will be started when it's its turn. If the service is under heavy load, then its turn might come later than you expect. Then, if the service is under heavy load, even after launching your job, it might take longer to finish because it won't have enough resources to run quickly.
Since you say that the same job finished quickly on a new tab, it sounds like either the first job was last in the queue or there was just a bug that caused it not be launched or something like that. We can't really help with things that resolve themselves like that.


text:
TBLASTN error: Too many positional arguments
I was running tblastn using standalone blast2.9.0+ through Ubuntu Linux, to get the alignment result of a list of protein sequences for same species, against its complete genome file.
The command I was using:
~/softwares/ncbi-blast-2.9.0+/bin/tblastn -query TBLASTNstring/882.Tstring -subject ../com_gen/*882*.fasta -outfmt "7 acc qacc pident qcovs evalue sstart send sframe" -max_hsps 1 -num_alignments 1 -out 882.blastout
(Info update: The query file has more than 1 sequence but all are in fasta format; the subject file has only 1 sequence; there is only 1 file with "882" named in the folder)
And I received error from tblastn as follow:
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: ../com_gen/DP_882_DeVul_Hilde.fasta

All the sequences in both my query and subject sequence files are of two lines, the header line started with ">", followed with sequences.
I have checked both query and subject sequence files using sed -n '0~2p' | grep -v -E '[[:alpha:]]' to see if there exist special character or space within the sequence line, but the output of the results return nothing.
PS: the complete genome file contained of nucleotide sequences with only ATGC (I have checked using grep -v "ATGC")
I ran the same command for other species as well and it worked fine.

Info Update: The screenshot above was example of the result for species with TaxID 196627.
I wondered what is the thing that went wrong here? And what else I can check on them?

How many files match the ../com_gen/*882*.fasta glob? Can -subject take more than one file as an argument? Also, aren't you mixing query and subject? What is TBLASTNstring/882.Tstring? That doesn't look like a fasta file, could that be the target DB name?

-subject and -query only allow one fasta file as argument

@PeterMenzel I haven't run blasts in ten years now. Do the modern versions really accept fasta files as -subject? Do they create the DB on the fly? Don't you need to first build the DB from the fasta file?

Hi @terdon I'm sorry for not being clear enough in the description. I have just updated some info in the post which is relevant regarding your question. There is only one sequence file named ../com_gen/*882*.fasta in my folder. (I was using * as the file name is quite long but certainly there is only 1 sequence file with the name of "882".) TBLASTNstring/882.Tstring is my query file which has more than one sequence but all in fasta format.

@PeterMenzel As far as I am concerned, blast can detect all sequences in a file given if they are in fasta format? Please correct me if I am wrong.

Anyway I just tried extracting one of the fasta sequences from the query file and ran it again against the genome. The same error still appear.

source: Stack Exchange

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

sort: votes


text:
Single or partitioned model for supermatrix tree inference?
I am dealing with several genes in my dataset. For each of the genes, I have built gene tree, with their best estimated model. I am intended to look into the effect of the tree from concatenated sequences (supermatrix tree), using Maximum Likelihood approach.
I am considering two options to apply the model for the supermatrix tree inference.

Obtain a single best estimated model for the supermatrix gene
Apply partitioned model for the supermatrix gene, based on each gene's model

From what I have read in most if not all of the research articles, they seemed to be using option 1, where a single model applied for the supermatrix tree inference. Since the purpose of using evolutionary model is to describe the changes in the gene and given that the evolution pattern in each gene differ, why wouldn't the respective model of each gene be considered in the tree inference?
For examples,
Dataset A: gene1 and gene2 are having LG and WAG model respectively as the best estimated model. But when the two genes are being concatenated, it turned out that LG is the best model for the matrix. It seemed that most if not all of the experiment opts for option 1 in this case. Theoretically, will the supermatrix tree using option 1 and option 2 in the model application still similar? Also, will the evolutionary pattern in gene2 still being represented correctly on the tree when we opt for option 1 in the supermatrix? Or it actually does not impact much on the calculation of the tree?
Dataset B: gene1 and gene2 are having LG as their best estimated model. The concatenated matrix has LG as its best estimated model too. Will model partitioning has any effect in this case? Does the ML tree/site likelihood calculation affected by the length of the input sequences?
Updated (replied from the RAxML-NG author):

combining multiple genes with the same substitution model into a
single partition would result in joint estimation of free model
parameters (e.g., alpha parameter for +G). So in general, it will make
a difference for likelihood computation, except if there are no free
parameters to estimate.

I understand concatenation well. What I'd like to be clear about is the 'super-matrix', is this the single mutation matrix derived from the concatenated alignment? Secondly are you obtaining your matrix before you perform the tree calculation and then inputing this matrix for the ML tree search?

@M__ yes for both. Say if I have 5 genes. The supermatrix derived from adding gene 2 at the end of gene 1, gene 3 at the end of gene 2, ... , gene 5 at the end of gene 4. They will exist as a big and single concatenated alignment matrix, concatenating gene1 to gene5 side by side. Then the supermatrix would be used together with the model/s to calculate for the ML tree.

Thanks, two final questions are these genes physically linked? Or random amplicons, across a genome? Finally, is this a mix of rRNA genes and protein genes and within the protein genes are they all single copy house keeping genes? I'll answer before the end of the week.

The genes are random across genome and yes they are  all single copy rRNAs and protein genes.

Thanks. Sorry to be persistent, are the protein genes 'house-keeping', e.g. often a core metabolic enzyme and all part of the core-genome, specific to the Order/family/species or a mix of the two?

No worries. It's just that my question is actually referring to those in general situation. It would be great if you could help me understand the concern. Probably I could justify it clearer.

Okay, you can't concatenate genes into a singular matrix particularly if these are rRNA genes versus protein genes, except in rare situations where the tests say 'thats okay' (below). This is simply because they evolve at very different rates rRNA genes have loci which evolve very slowly and are highly conserved.
The theory goes back to incongruence length difference test (ILD), which is somewhat confusingly also known as the partition homogeneity test and was then modified to become the localised ILD test.
The modern version of this test, simply called 'homogeneity' is found in IQ-Tree. Thus, the best approach is to test the assumption empirically,

IQ-TREE provides three matched-pairs tests of symmetry (Naser-Khdour
et al., 2019) to test the two assumptions of stationarity and
homogeneity. A simple analysis:
iqtree2 -s example.phy -p example.nex --symtest-only

Just to be clear this tests whether the rate between two partitions (genes) is homogeneous. If 'yes' the supermatrix is fine, if 'no' the data needs partitioning.
The theory is simply that if two genes evolve at noticeably different rates they can't be combined into a singular matrix. The theory goes back a long way is used to circumvent long-branch attraction issues. Thus, at a minimum the rRNA genes and protein genes need to be compartmentalised, i.e. have their own matrix. Between protein genes, if they are house-keeping genes, generally they'll have homogeneous rates (which you can empirically test via IQ-TREE).
If someone is concatenating rRNA with protein genes without accommodating rate variation via partitioning it is simply sloppy phylogenetics and if published represents a shortcoming with the reviewing process. This is because in all certainly these loci will fail any rate homogeneity test, there's a slim off-chance its okay which you can test it empirically. Say if 18S and 28S were concatenated (or just 16S if these are bacteria) with an atp gene, that can give a different tree to 18S with atpA, atpB, atpC, atpD etc... so the ratio between rRNA to protein genes can start influencing the tree ... which is simply wrong.
If someone concatenates a load of house-keeping genes thats okay because there's a good chance they're homogeneous and if even the test failed, they may not be that far out, rRNA versus proteins can be way out. This can easily be an order of magnitude out and thats difficult to ignore. When there's proteins under selection versus conserved housekeeping genes ... its very unlikely homogeneity would hold.
Summary As a general rule of thumb in your method of tree building compartmentalising all genes is safer route, computationally there's no difference between zero partitions and a partition for every gene. You are never going to make a mistake.
Where it gets complicated is when the phylogenetic parameters, particularly in advanced models, are themselves estimated via maximum likelihood alongside the topology. The rational for this is complicated. Computationally, this can get very, very expensive even on supercomputers and hence there's a tendency to make an assumption (often without checking it) that there are fewer partitions than there actually are. Regardless, assuming away rRNA versus protein genes even in this method of tree estimation should not done and the test will almost certainly fail.

Appreciate your answer. I agree with the possible induced error when concatenating together rRNAs and protein genes. Sorry for the confusion caused from the comments. What I mean was I have dealt with both type of datasets and they were concatenated separately. But that was quite an interesting input for me, especially the ILD test, will try that out.

Btw, from your last paragraph, though not really familiar with the concept behind, I was aware that it is possible to estimate best models via the ML algorithm. But what does it meant by "there's tendency to make an assumption that there are fewer partitions than there actually are"? Also, I am very interested with "computationally there's no difference between zero partitions and a partition for every gene", would you mind to elaborate a bit more if possible? What do you think about it theoretically? Btw, I have updated the question. Thanks!

Agreed @web the error was unfortunate. "computationally no difference ...", an exaggeration, should be "the difference in computation time is not a major issue". Final point "tendency ...", a reviewer will be aware of the difference between what is computationally practical against the theoretical requirements of the calculation, certain approximations are considered ok.

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


text:
Replicating VCF Filtering & Trait-Based SNP Extraction Workflow
I'm new to handling SNP genotyping data in VCF format.
My goal is to identify significant SNPs linked to a specific trait (like "height") for specific samples in a multi-sample VCF file. My previous colleague who initially handled this project left, and I'm struggling to replicate their workflow for filtering and trait-based SNP extraction.
Here's what I'm working with: a large VCF file containing data from 32 samples, ideally needing a filtered subset focused on just 5 specific samples. These VCF files are GATK outputs with all detected SNPs. I want to extract statistically significant SNPs that are non-homogeneous (meaning they vary) for the target trait (height).
GATK command in the VCF file (If this helps):
##GATKCommandLine=<ID=SelectVariants,Version=3.3-0-g37228af,Date="Thu Jan 12 15:24:24 CET 2017",Epoch=1484231064776,CommandLineOptions="analysis_type=SelectVariants input_file=[] showFullBamList=false read_buffer_size=null phone_home=AWS gatk_key=null tag=NA read_filter=[] intervals=null excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=/work/NGSwaiting4newNAS/G01_chromosome.fna nonDeterministicRandomSeed=false disableDithering=false maxRuntime=-1 maxRuntimeUnits=MINUTES downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=1000 baq=OFF baqGapOpenPenalty=40.0 refactor_NDN_cigar_string=false fix_misencoded_quality_scores=false allow_potentially_misencoded_quality_scores=false useOriginalQualities=false defaultBaseQualities=-1 performanceLog=null BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 globalQScorePrior=-1.0 validation_strictness=SILENT remove_program_records=false keep_program_records=false sample_rename_mapping_file=null unsafe=null disable_auto_index_creation_and_locking_when_reading_rods=false no_cmdline_in_header=false sites_only=false never_trim_vcf_format_field=false bcf=false bam_compression=null simplifyBAM=false disable_bam_indexing=false generate_md5=false num_threads=1 num_cpu_threads_per_data_thread=1 num_io_threads=0 monitorThreadEfficiency=false num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false variant_index_type=DYNAMIC_SEEK variant_index_parameter=-1 logging_level=INFO log_to_file=null help=false version=false variant=(RodBinding name=variant source=DBG_SNP_EGRESQ_filtered.vcf) discordance=(RodBinding name= source=UNBOUND) concordance=(RodBinding name= source=UNBOUND) out=org.broadinstitute.gatk.engine.io.stubs.VariantContextWriterStub sample_name=[] sample_expressions=null sample_file=null exclude_sample_name=[] exclude_sample_file=[] select_expressions=[] excludeNonVariants=false excludeFiltered=false restrictAllelesTo=BIALLELIC keepOriginalAC=false mendelianViolation=false mendelianViolationQualThreshold=0.0 select_random_fraction=0.0 remove_fraction_genotypes=0.0 selectTypeToInclude=[SNP] keepIDs=null fullyDecode=false forceGenotypesDecode=false justRead=false maxIndelSize=2147483647 ALLOW_NONOVERLAPPING_COMMAND_LINE_SAMPLES=false filter_reads_with_N_cigar=false filter_mismatching_base_and_quals=false filter_bases_not_stored=false">

Further complicating things, I found previously filtered VCFs with significant SNPs for 5 samples. Now, I'm tasked with repeating the process for more samples.
While I've researched VCF filtering by origin selection, I'm stuck on pinpointing the most efficient workflow to achieve both filtering and trait-based SNP extraction. Any recommendations for tools or resources would be incredibly helpful! Please feel free to ask if there is any elaboration needed. Thank you!
Update: I also have another file with a list of QTL regions that were significant to the height trait.
Below is a picture generated to help with understanding of the challenges:

I have information on the specific region of QTL that are significant to the height trait. But the information is based on all the available origins. Now I am only interested to some of the origins. How do I filter the VCF files to get the SNP related to the interested origins?

What does "origin" mean? Do you mean different individuals were sequenced and all results are in a single, multi-sample VCF? And do you only have SNPs? Nothing larger? Do you have a budget or do you need to only use free tools?

@terdon Yes, origin refers to different sample included in the genotyping and all the results were compiled in one VCF. I also have the QTL regions that are significant to the trait, the fasta sequences for each chromosomes, etc. I think the files are complete as they have completed the task previously. I would like to explore the free tools if possible.

OK, and what do you have in the VCF file? Is that annotated with the gene the variant falls on/near? Do you have a list of genes for your phenotype? Again, are these only SNPs (normal GATK output would also include indels)?

extract statistically significant SNPs that are non-homogeneous (meaning they vary) for the target trait (height).  assuming controls (A, B, C) and samples with the desired trait (X, Y, Z) what does it mean? In the order of the letters above, anything different from say 6 Ts, ABC=reference base, XYZ=T, ABC=whatever, XYZ=CGT?

@terdon The SNP is not annotated with the gene, I think the information is available somewhere else, but not in the VCF file. At the moment, I can only see list of QTL significant to the phenotype and several sequences related to the phenotype from publication. I am not sure how to check if indels were included in the VCF file, is it like what explained by Martijin here (https://www.biostars.org/p/48731/#48736)? If yes, then no, as the ref and alt allele consisted only a single base pair here?

@darked89 From what I understand, in the VCF file, there are a reference and alternate allele. The genotyping result for each sample are based on this. For example the genotype received by sample A is 0/0, sample B is 0/1 and sample C is 1/1. In this case, sample A is having G/G, sample C is having T/T, which are homogeneous, only sample B is having G/T, the non-homogeneous I mean in the text.

@terdon are the indels usually being noted in the "ID" column? My "ID" columns are filled with only dot. I have updated the question with GATK command I found in the VCF file.

No, the ID column will be empty unless you've annotated the VCF with a tool that populates them. Indels appear as you said: >1 residue in either the REF or ALT fields, or both.

I think what you are looking for is to perform a GWAS: Genome-Wide Association Study. Generally, these are used to detect links between mutations in patients and complex diseases. I see no reason that the same approach couldn't be used to attempt to predict someone's height.
I haven't run a GWAS in a long time, so I wouldn't be able to recommend any specific tools to use, but I found packages in R and in Python that could work.
One thing you didn't mention is how you know the height of each of your donors, which would seem like the biggest complexity to me. Once you have the group VCF + associated NGS files, you're mostly covered from a data standpoint.
If you're looking to go deeper, you could investigate dimension reduction methods like PCA that can also serve to associate variables with certain traits, or even some machine learning methods like SVMs or a Neural-Net regression model.

Thanks for the response! My target is actually plant, it was the plant height I have been referring to, but I guess that is not making much differences? There were several samples involved in the sampling and genotyping. The height of plants were measured at nth years. I am confused with what you have mentioned: "group VCF + associated NGS files", could you elaborate more on this?

GWAS should work on plants just fine, even if it was designed for use on Human diseases (these two articles use GWAS on plants https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8316211/ https://www.nature.com/articles/s41396-021-00993-z maybe their methods detail their parametres). The "group VCF" refers to the file containing the variant calls for your entire population, which can be quite long to obtain and requires a lot of computing power, so having access to it is a boon. If you also have access to the raw NGS data, then you can also explore other methods of analysis, like k-mers, if you want.


text:
Applying glmnet to identify predictors for subtypes
I am new to glmnet and other ML techniques, so apologise in advance if this sounds a trivial question. However, your guidance is very much appreciated.
So, I have nearly 700 genes with their expression profile across 4 conditions. I want to identify signature genes for each of the 4 conditions. I am using glmnet but I am not getting anywhere near the results.
so my data looks like this
           Subtype     HES4   FNDC10    PRKCZ PRKCZ.AS1
Sample1     NRT 5.452955 6.241387 8.546059  6.260245
Sample2     RT 5.125350 7.038719 8.900024  6.267325
Sample3     NRT 5.689560 5.903811 8.029434  4.213996
Sample4     NRP 5.834468 5.722589 8.061339  5.072274
Sample5     RP 5.401623 7.272111 9.108342  6.599629

purely in terms of numbers, the subtypes have following number of samples.
NRT 52
RT 42
NRP 30
RP 7

Q1. May be the groups are imbalanced and caution should be taken?
I have applied glmnet as following
modellinglasso=read.table("LRMatrix.txt", header=T, row.names = 1, sep="\t" )
modellinglasso$Subtype=as.factor(modellinglasso$Subtype)
modellinglasso$Subtype
#build a single pass (single fold) lasso-penalised model
lassoModel <- glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  family="multinomial")
plot(lassoModel, xvar="lambda")
Warning message:
In lognet(xd, is.sparse, ix, jx, y, weights, offset, alpha, nobs,  :
  one multinomial or binomial class has fewer than 8  observations; dangerous ground

#perform 10-fold cross validation
cv.lassoModel <- cv.glmnet(
  x=data.matrix(modellinglasso[,-1]),
  y=modellinglasso$Subtype,
  standardize=TRUE,
  alpha=1.0,
  nfolds=10,
  family="multinomial",
  parallel=TRUE)
# plot variable deviances vs. shrinkage parameter, λ (lambda)
plot(cv.lassoModel)

idealLambda <- cv.lassoModel$lambda.min
idealLambda1se <- cv.lassoModel$lambda.1se
print(idealLambda); print(idealLambda1se)

[1] 0.1186125
[1] 0.1186125

# derive coefficients for each gene
co <- coef(cv.lassoModel, s=idealLambda, exact=TRUE)
head(co)

co.se <- coef(cv.lassoModel, s=idealLambda1se, exact=TRUE)
head(co.se)

# identify predictors for each sub-type
rownames(co.se$PR)[which(co.se$NRT != 0)]

[1] "(Intercept)"

output for
co.se$RT!=0
6 x 1 sparse Matrix of class "lgCMatrix"
            1
(Intercept) |
HES4        .
FNDC10      .
PRKCZ       .
PRKCZ.AS1   .
MMEL1       .

It doesn't matter which subtype I put, I get the same output "(Intercept)" instead of gene names.
Any help to get the gene signatures for the subtypes is sincerely appreciated.
Thank  you

Final answer. Keep in mind I work in Python so I'm trying to translate here.

The output of coef(cv.lassoModel) ... those are your genes of interest thats your feature selection. head I assume gives you the top values ... I think thats right (I need to check its not coef(lassoModel)) ... Basically thats where your results are.

This might be over simplistic but co.se without the head is your answer. Not every gene will carry a weight BTW. If the entirety of that is (700) . ... try point 2 below and then point 4. I don't think that will happen (see Whats happened)

I would redo this with alpha=0 (ridge regression) and alpha=0.5 elasticnet, but you don't need to
Generate an accuracy score

    y_predicted <- predict(best_model, s = best_lambda, newx = x)
    sst <- sum((y - mean(y))^2)
    sse <- sum((y_predicted - y)^2)
    rsq <- 1 - sse/sst
    rsq

Hopefully this is large e.g. around 0.8 (in Python this automatic BTW). You could ignore that, but its what I'd do

The intercept is simply the weight when your gene expression score is 0.
Lasso regression handles sparse data - so no worries about ".". I don't whether that Python's version of "nan" (not a number), or if its zero.  ... @haci knows what this means and his answer is authoritative for the second part of your question.
In Python transformations are common, but that doesn't seem an option here, so find providing your accuracy score is reasonable, personally I'm happy.

In Python what you want is feature selection
print(pd.DataFrame({'Features': X_train.columns, 'Coefficients': lasso.coef_}))
Whats happened? The coefficients for head(co.se) have been shrunk to zero, thats only the top 10? results, you've 700 weights. The full list should (famous last words) contain a minority of really big coefficients ... those are really important and those are your answer. So by printing only the top few results ... not the whole lot the results are still there you've just not seen them. Thats it, basically you were almost there and you could happily settle on that as the final answer. Basically, alot of genes ain't doing much a few genes are doing all the work - thats the answer. Below I stress repeatedly I don't think anything is wrong here and I personally like the result.

Background
Lasso regression, it's a cousin of ridge regression and key parameter is lambda (who'd guessed). What your plots are doing is tuning the parameter and finding the associated weight individually or collectively. Fairly straight forward. So like whats the value of lambda for your data set ('answered' below)?  It's worth doing some reading on how lambda is calculated to understand the plot, if you're not sure how it works.
Imbalance, well to be honest I'm not convinced there's anything wrong with your analysis. Okay there's a lack of stationarity for some of your weights ... but some are fine and gist is pretty clear. Yeah okay there's a problem with stabilisation of lambda - thats not necessarily imbalance.
Imbalance is when it's 1 in 10 and over. You've got 1 in 6, 1 in 7, not huge. However, - you might be right - and if you're worried simply repeat the analysis excluding RP (value 7). Then check whether the parameterisation plot (lambda against weight) is more stable or more "homogeneous" (multinomial plot), i.e. the genes behaviour in a more uniform way. If that is does work - it gives you nicer results AND RP is essential to the analysis thats a separate question about augmentation. That is complicated. For example, the multinomial plot with a 3 target categorical model might give a clear value for lambda.
Transformation Here, I would suggest if you want to stabilise the multinomial plot (if it's possible) I would look at a transformation. The data you presented doesn't look like it needs a transformation, but we can't see the whole data. In Python there's a max-min transformation - thats where I might start.
Overall good results Personally, I'd say these were good results, except the weights ain't homogeneous for some genes (most are fine BTW), but the patterns clear. The red (top), blue (bottom) and cyan (negative y-axis) are real results showing that within the 700 genes 3 genes are having a massive impact. Sure you have positive and negative and thats how it works one gene is positively contributing to the classification, others are negatively contributing to it. Then you've the black and 'deep pink' lines which have stabilised around 5E4. From 700 genes down to 5 genes. Whats wrong with that?
Whats the answer to the data set?
Its ...

NRT
RT
NRP
RP

... is due to 'red', 'blue', 'cyan', 'black' and 'deep pink' genes. Thats it.
Frank assessment To be perfectly honest, looks good. No-one is going to say "Wooaaa lambda ain't stable, destroy! ;-)", they'll say "Wow thats cool and clear".

How I'd do it.. I would look at basic ML robustness statistics and use this to assess the model rather than the parameterisation of lambda. Its horses for courses as they say. Upvotes.
Note I use Python not R for ML, both work well however. glmnet within ML is new to me BTW.

Thank you for your reply and detailed explanation. What I do not understand is why I am getting output as Intercept instead of gene names.

Another puzzling aspect is when I do co.se$RT!=0, I get a gene list with "." next to it :(

Technically you are getting a sparse Matrix with co.se$RT!=0 and the "dots" are representing zeros or in this specific case FALSE values. You can see yourself by as.matrix(co.se$RT!=0).

All of them :( I tried to normalise the values by getting zscores. HES4 is a gene name and its reading in different samples.

@Angelo I've addressed the comments in the first part of the answer.


text:
shiny code with seurat object
I am trying to create a shiny board for my scRNA seq data. I have saved my seurat objects adn I use following code as below:
library(shiny)
library(shinythemes)
library(Seurat)
library(uwot)

set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")

ui <- fluidPage(

    # Application title
    titlePanel("scRNAseq data browser"),

    fluidRow(
        column(4,

               sidebarPanel(
                   radioButtons("dataset", label = h3("Dataset"),
                                choices = list("sample1" = "set_1", "sample2" = "set_2"),
                                selected = "set_1"),
                   helpText("Enter gene names in CAPS.  Gene names must be exact."),
                   textInput("gene2", label = "Gene Name", value = "GAPDH")

               )
        ),
        column(8,
               column(6,
                      mainPanel(
                          plotOutput("dimPlot1", width = "150%"),
                          plotOutput("genePlot2", width = "150%")
                      )
               )
        )
    )
)

server <- function(input, output, session) {

    datasetInput <- reactive({
        df <- (input$dataset)
            })

    output$dimPlot1 <- renderPlot({
        DimPlot(datasetInput(), reduction = "umap")
    })

    output$genePlot2 <- renderPlot({
        FeaturePlot(datasetInput(), features = (input$gene2), reduction = "umap")
    })

}

shinyApp(ui, server)

That error that I get on the section where image needs to be displayed is "Subscript out of bound".
Following dataset can be downloaded and used as an example curl -Lo pbmc_10k_v3.rds https://www.dropbox.com/s/3f3p5nxrn5b3y4y/pbmc_10k_v3.rds?dl=1
Can anyone help in fixing this shiny script so that it works.
Many thanks.

Not easy to reproduce your error but here is my guess. This is a string, you can check it:
input$dataset
class(input$dataset)

This is a simple solution. If you want to call out the dataset depending on input, use a list:
set_1 <- readRDS("/temp2/data/mydataB10.rds")
set_2 <- readRDS("/temp2/data/mydataB5.rds")
data_list = list(set_1=set_1,set_2=set_2)

Then in the server :
datasetInput <- reactive({
        df <- data_list[[input$dataset]]
            })

It worked, thank you:)

@StupidWolf if you can give a look at this https://bioinformatics.stackexchange.com/questions/21043/shiny-app-for-single-cell-cluster-export-issue


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

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

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

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

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

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

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

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

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

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


text:
What pitfalls exist with running FastQC on a bam file?
I was analysing a bad sequencing run of some RNA data using FastQC, I supplied it RNA-STAR-aligned bam files (and the human reference).  The output of MultiQC had many more "unique reads" in the sequence counts table than there were reads in the FastQ files that were fed to RNA-STAR.
At first, I realized that all the FastQC stats were about "mappings", not reads (as they were labeled), but "unique mappings" doesn't seem to make sense either, because even if I consider the possibility that segments of reads were mapping in different places, the numbers don't make sense.  Here is the worst case:

Max read length after trimming: 47 bases
Number of reads in the original FastQ file that went into STAR: 59M

Bam Stat Table Numbers:

Total records: 391 million
Unique mappings: 3.1 million
Non-primary Hits: 335 million

FastQC Sequence Counts Bar graph:

"Unique reads": 185 million
"Duplicate reads": 206 million

While I found documentation stating that you can run FastQC on bam files, there were no caveats I could find about it needing to be unmapped reads.
Can anyone explain precisely why there can be 185 million "unique reads" when bam stats says there were 3 million unique mappings and the number of reads is 59M?

I don't see why you are trying to compare two totally different measurements.
When FASTQC says a read is unique, it means that no other read shares its sequence.
When STAR says a read is uniquely mapped, that means it maps to one and only place.  That does not mean that there aren't a hundred more identical reads that map nicely to the exact same place.
You could have a very high level of PCR duplication of reads that nicely align to one and only one place. Or you could have a whole bunch of reads that differ slightly from each other who all map to 5 locations each.
FASTQC was designed for DNA sequencing, not RNASeq.  Some of its warnings make no sense in the context of RNASeq.  If I was doing RNASeq of a filtered set of pancreatic alpha cells,  I would expect a hell of a lot of reads for glucagon, and a lot of them would be natural duplicates, and that would be biologically normal, and not a problem.
You could take the same batch of reads, and align them to one genome, and align them to a mix of two closely related genomes.  The second will have far more non-uniquely mapped reads, but the FASTQC stats will not have changed.
That said, if you doing human or mouse RNASeq, 1% unique alignment means something is wrong.  Find out where all those reads are aligning, you might have rRNA contamination.  If you are doing an organism with a high level of duplication in its genome that might be normal.

But the number of reads in the FastQ sample file that was fed into STAR was 59,194,421 (I stated it was much fewer than was reported).  I believe that the 185M and 206M (391M total) referred to by FastQC as "reads" are actually "mappings" since I fed FastQC a bam file instead of a FastQ file, which is what this question is about. Perhaps this is a misunderstanding stemming from the recent edit. I will fix it.

I understand the general caveats about running FastQC on RNA-Seq FastQ data. I regularly check for rRNA contamination and check over-represented sequences for rRNA segments. My question however pertains to why you can run FastQC on a bam file, what its utility is, how to interpret the results in that context (e.g. a "read" is a "mapping"), and what the pitfalls are WRT using a bam as input. Also, I knew it was a bad sequencing run, but this question is not about identifying bad sequencing runs.

Hi I think @swbarnes has answered the question very nicely, they are concerned with two different measurements and swb proivide the underlying rationale ... one was for DNA genomes not for RNAseq. It would be cool  to accept the answer.

If you think that FASTQC should have a built-in filter for primary alignments only when reading a bam, tell the FASTQC developers, not me.


text:
What column and row naming requirements exist with Seurat (context: when loading SPLiT-Seq data)
I'm trying to use Seurat for the first time and am learning about single-cell analysis for the first time, and I'm doing so with split-seq data. (Full disclosure: I'm also a lightweight when it comes to R.)
When Read10X didn't work on the mtx and associated files output by the split-seq-pipeline no matter how much I tried to make it look exactly like a 10x dataset (based on numerous examples and descriptions), I tried converting my sparse market matrix format into a dense tsv format (see update below for a sample of the file) and I read it in using:
dge_data <- as.matrix(read.table("DGE.tsv", sep="\t", header=TRUE))
pbmc <- CreateSeuratObject(counts = dge_data, min.cells = 3, min.features  = 200, project = "SPLITSEQ_TT", assay = "RNA")

I started following chapter 9 of the scRNA-Seq course I found here. I ran into some issues which I was able to resolve by tweaking my mtx to tsv conversion script, such as including the chromosome in the row names (so that the mitochondrial QC step would work), but I ran into some errors that I assume have to do with the format of the gene names/IDs and possibly some issues stemming from cell names/IDs (which I've thus far been able to work through I think).
However, given errors like the following, and the fact that I can't seem to find specifications on these sorts of row/column name requirements (e.g. no commas and must include at least the "MT" chromosome name for QC steps using mitochondrial genes), I was wondering if I could be directed to somewhere where those specifications are laid bare. For the most part, the gene names have not been a problem, but I run afoul in this step (where it appears to complain about a comma):
> FeatureScatter(object = pbmc, feature1 = "hg38-4,ENSG00000074966,TXK", feature2 = "PC_1")
Error in parse(text = x) : <text>:1:7: unexpected ','
1: hg38.4,
          ^

Also, when I look at how CellsByIdentities behaves (while I was working through other errors/warnings), I realized I'm not sure that the column names are being processed correctly due to the underscores (followed by a number (e.g. _0)), because my understanding is that the _# portion of the split-seq cell barcode indicates the first well barcode, which would be needed to uniquely identify the cells in the table. I had retrieved the cell barcode IDs from the cell_metadata.csv file output by split-seq. Should I modify these column names? Because I infer from this output from CellsByIdentities that it's only using the sequence portion of the cell barcodes to uniquely identify cells, which I don't think is correct:
> CellsByIdentities(object=pbmc, cells=colnames(pbmc))
$AAACATCGAAACATCG
[1] "AAACATCGAAACATCG_0" "AAACATCGAAACATCG_1"

$AAACATCGAACGTGAT
[1] "AAACATCGAACGTGAT_0"

$AACGTGATAAACATCG
[1] "AACGTGATAAACATCG_0" "AACGTGATAAACATCG_1"

$AACGTGATAACGTGAT

[1] "AACGTGATAACGTGAT_0" "AACGTGATAACGTGAT_1"

The other assumptions have to do with data size. I was running a small set of test data, and I had to adjust some parameters. So I would like to know if there's a hard lower limit on the data. And do I need to restructure my data to separate counts of genes from different species (e.g. mouse & human)? I.e. Should I create separate matrices that have only the human-identified cells and gene annotations to use it with Seurat?
I can't seem to find these specifications, presumably because Seurat was designed exclusively around 10x data and that the tool wasn't tested on data from other sources.  Perhaps I'm worrying too much, but without knowing what may rely on column/row name format assumptions, I can't be sure my data is being analyzed properly.
UPDATE
Here is the head of the tsv file I originally tried to load when I ran into various errors that I suspect were due to the formatting of the row and column names (note, there are only 7 cells total, and it's a small but real test dataset from our core facility):
AAACATCGAAACATCG_0  AAACATCGAAACATCG_1  AAACATCGAACGTGAT_0  AACGTGATAAACATCG_0  AACGTGATAAACATCG_1  AACGTGATAACGTGAT_0  AACGTGATAACGTGAT_1
hg38_X,ENSG00000000003,TSPAN6   8   2   2   26  15  17  13
hg38_X,ENSG00000000005,TNMD 0   0   0   0   0   2   0
hg38_20,ENSG00000000419,DPM1    22  26  27  73  43  50  28
hg38_1,ENSG00000000457,SCYL3    6   6   5   12  7   19  8
hg38_1,ENSG00000000460,C1orf112 23  10  24  59  34  49  21
hg38_1,ENSG00000000971,CFH  101 44  129 331 151 270 122
hg38_6,ENSG00000001036,FUCA2    33  28  24  82  58  74  36
hg38_6,ENSG00000001084,GCLC 30  16  14  72  27  54  42
hg38_6,ENSG00000001167,NFYA 14  9   17  41  28  31  12

Yesterday, I did end up trying to change the column and row names to the following:
AAACATCGAAACATCG-0  AAACATCGAAACATCG-1  AAACATCGAACGTGAT-0  AACGTGATAAACATCG-0  AACGTGATAAACATCG-1  AACGTGATAACGTGAT-0  AACGTGATAACGTGAT-1
hg38-X.ENSG00000000003.TSPAN6   8   2   2   26  15  17  13
hg38-X.ENSG00000000005.TNMD 0   0   0   0   0   2   0
hg38-20.ENSG00000000419.DPM1    22  26  27  73  43  50  28
hg38-1.ENSG00000000457.SCYL3    6   6   5   12  7   19  8
hg38-1.ENSG00000000460.C1orf112 23  10  24  59  34  49  21
hg38-1.ENSG00000000971.CFH  101 44  129 331 151 270 122
hg38-6.ENSG00000001036.FUCA2    33  28  24  82  58  74  36
hg38-6.ENSG00000001084.GCLC 30  16  14  72  27  54  42
hg38-6.ENSG00000001167.NFYA 14  9   17  41  28  31  12

And it seems to have fixed most of the issues.  The violin plots now have curvy brown swaths down the center that were not there before.  And in particular, when I run CellsByIdentities manually on the column names now, I get what looks correct (as opposed to what I pasted yesterday above (which appeared to group the 7 cells in 3 pairs and one individual):
CellsByIdentities(object = pbmc, cells = colnames(pbmc))
$SPLITSEQ_TT
[1] "AAACATCGAAACATCG.0" "AAACATCGAAACATCG.1" "AAACATCGAACGTGAT.0" "AACGTGATAAACATCG.0" "AACGTGATAAACATCG.1" "AACGTGATAACGTGAT.0"
[7] "AACGTGATAACGTGAT.1"

I think that the remaining issues are likely due to the small data size.  For reference as to the data size, this is the first data line from the original mtx file:
7 38526 194539

Please note also that I'm not looking for how to get this particular dataset to work.  I am (as a member of the bioinformatics team that supports users of our core sequencing facility) trying to establish the requirements for the conversion script so that when we deliver SPLiT-Seq results, our output can be directly used with Seurat by the users without issues or need for further conversion.  So I specifically need to know what makes row or column names valid or invalid, as opposed to what to change WRT this specific dataset.  We may be supplied with annotations from various sources that may have different characters in the gene symbols, or accession numbers, or even chromosome names.

hepcat72 you have thanks @haci for their answer. Please do consider "accepting" the answer, because it appears close to what you are requesting and/or upvoting.

Read10X() works on mtx files, that is not the root of your problem. The problem is that you are pushing non-10x-data to look and work like 10x-data.
You can just use CreateSeuratObject() and feed your count matrix to this function. I am not sure if this particular function accepts mtx files, could not find this in the documentation, but in the worst case you will just need to supply a plain count matrix, count matrix being genes as rows and cells as columns.
And that could be a good idea to show the first few rows and columns of your SPLIT-Seq output so that we can have and idea.
Regarding gene identifiers, I would stick with gene symbols, to help you in cluster annotation but also using the "MT" prefix for mitochondrial genes for example.
It looks like you are providing somewhat obscure identifiers as gene identifiers, a mix of genome build, ensemble gene id and gene symbol and apparently Seurat does not appreciate the ,. While I cannot give you a definitive recipe on what can and cannot be set as gene identifiers when using Seurat, I can tell that I have worked both with gene symbols and Ensembl identifiers, separately, and both was fine with Seurat.

Yes, I converted my data to tsv where rows are genes and columns are cells. Setting aside my concerns about cell column names, don’t the row names have to be unique, especially when the genes are a collection from 2 species? I also noted that the hg38 Ensembl gene symbols already have MT- in them, but mm10 doesn’t, so I would have to address that.

In theory the human and mouse gene symbols would not be the same for the "same" (actually orthologous genes); human gene symbols would be all uppercase and those of mouse would start with an uppercase letter followed by lowercase. Anyway, you can alter gene names anyway you like, doing that before reading the data into Seurat would be the "safest" I believe.

And yes, I’m feeding CreateSeuratObject the dge_data I described above.

I also did end up modifying the column names (the cell barcode IDs) to replace the underscore with a dash. When I do that and call CellsByIdentities, it returns a single list of 7 cell barcodes instead of 3 pairs and 1 single barcode, so it seems like not having underscores in those column names is critical? I tried to address all of the row/column name issues I encountered and after doing so, the plots all look much better. For example, the violin plots now all have curvy brown columns that weren’t there before I changed the row/col names for commas and underscores.

I have edited my answer, basically I would not have commas in my gene identifiers. Secondly, I have underscores in my cell names and Seurat is fine with it, in fact I believe it was one of the Seurat functions appending underscores to cell names while merging data.

Thanks so much for your help. I will add a sample of the data to the question in a bit.

OK, I updated the question to add the requested sample.  Note, the core facility requested the cell barcodes to be used as the column names. Also, WRT hg38's and mm10's gene symbols' uniqueness... My background is more prokaryotic. Gene symbols in that realm are definitely not unique and my tentative understanding was that symbols for eukaryotes were not unique either, particularly when it comes to paralogues. However, I suspect that only unique mappings make it through the process, which mitigates uniqueness issues, but may not eliminate it totally. Does Seurat group by common gene name?

With the recent additions to your question, now it looks like it would be the best to address to the Seurat developers. They are quite active on GitHub.

Thanks. Yes, I posted a couple days ago (https://github.com/satijalab/seurat/issues/2887).  I guess I am just anxious to hear back. Your feedback has been valuable though because I wasn't sure whether I was doing something wrong or not. And it sounds like I am doing things roughly correctly and that these issues may not necessarily be due to me doing something wrong (though I'm still not entirely sure).


text:
Why doesn't an (Entrez eutils) einfo request for "gene" return the link gene_nucleotide or gene_nucleotide_pos links?
I'm updating a galaxy tool wrapper for Entrez's eutils suite and I'm trying to create a form with valid link selections (among other things) based on the "from" & "to" databases to reduce the size of the select list.
Einfo doesn't return these links for the gene or nuccore databases. I'd found this resource which seemed to show all the possible links, but it it doesn't show these links either. Yet, in working on elink's acheck feature, using some gene IDs in a test run, I saw 2 links that were not listed in either resource:

gene_nucleotide
gene_nucleotide_pos

I manually checked that I can use those as link names in valid elink queries. It did return valid results (i.e. no error).
So I ran elink on the gene database using a few random UIDs to see if I get back a comprehensive list of links (albeit possibly empty results for those links), but it also was missing these 2 links.  I'm concerned that I'm not getting all of the possible links for my form interface and that it could produce valid results and that users may want to link via these links, given prior knowledge that they exist.
How do I get a comprehensive list of database links from Entrez and why doesn't einfo return all possible links?

I've been corresponding with NLM about this issue and I finally took the time to try out their suggestion (which personally I found hard to see between the lines and which is not a discrete solution, but rather a very time-consuming manual process containing false positives because they say that to get a formal and comprehensive response to my query:

Collecting the linkname will be a difficult task, that will take time and coordinate/check with relevant groups maintaining individual databases. Your patience will be greatly apprecired. [sic]

It took me a bit to decipher their suggestion on how to manually find all links, so I will share what I've learned.  There is (currently) no codified or formal means to obtain all possible link names and I infer that many links exist simply because they are utilized in the internal workings of their website.  I could be wrong about that, but suffice it to say that there are many undocumented links and they change constantly.
You can get a list of filter items using each individual database's advanced search web interface, "most of which should be the linkname from that source database to other target database".  (So this should result in all possible links among a series of false positives.)
Here's how you do that.  Let's take the gene database as an example:

Go to the advanced search web page for that database, e.g. https://www.ncbi.nlm.nih.gov/gene/advanced
In the "Builder" section, select "Filter"
Click "Show Index List" to the right of the search term field next to the select list where you selected "Filter" (you will see a multi-select list appear below the search field)
Click "Next 200" until you find a selection in that new select list matching the database name followed by "all" and a number in parentheses, e.g. gene all (29994947) and click to select that item
Click "Refresh index" at the lower right corner next to the multi-select list

The multi-select list will be repopulated and most of the items will be linknames (though you have to replace spaces with underscores in order for those linknames to work with the elink utility).
Doing this for the gene database and then scrolling through just links to nuccore (and its alias "nucleotide" - ignoring est and gss), you will find:

gene nuccore
gene nuccore clust
gene nuccore mgc
gene nuccore pos
gene nuccore refseqgene
gene nuccore refseqrna
gene nucleotide
gene nucleotide clust
gene nucleotide mgc
gene nucleotide mgc url
gene nucleotide pos

If you look only at the linknames in the documentation or via einfo (by supplying the gene database and looking at similar links to nuccore), you only get these links:

gene_nuccore
gene_nuccore_mgc
gene_nuccore_refseqgene
gene_nuccore_refseqrna
gene_nucleotide_mgc_url

So I believe that given this information and the ever-changing nature of these links, I believe I will allow the user the option to enter a linkname manually, if the set in the select list I generate via einfo does not contain the link they need.

If you select the items in the resulting list, they are added to your search term in the above field where you can copy them out.


text:
Ribosomal RNA QC quantification using SILVA
I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.
I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually).  I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.
I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from SILVA would be great, but I don't know if I'm using the wrong fasta file or what.  I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.
So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found.  The others (5S, 5.8S, and 28S) don't appear to be there.  Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?
UPDATE: Please note, this question is about where to find rRNA and/or rDNA sequences to use as a reference, not about tools that can be employed in quantifying rRNA read content.  While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification.  Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool).  And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern.  It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.

Have you thought about using an existing tool to do this QC? RNAseQC (https://github.com/broadinstitute/rnaseqc) for example has has a line that tells you how many reads mapped to rRNA genes.

Using what reference & annotations? The 5S, 18S, and 28S are nowhere to be found. I couldn’t find them in any GTF files. I’m not an expert in rRNA, but I assumed they would be there. The researchers complained they weren’t, & when I looked, I was shocked I could not find them. I had to gather them up manually. I didn’t know about that output of rnaseqc, but I assume it’s dependent on the annotation you provide it.

I was using bwa in my galaxy workflow, but any relevant tool is fine. rnaseqc has a -BWArRNA <arg> which takes a fasta reference for rRNA sequences. It mentions that if that flag is absent, it uses the GTF file, which is says is fast, but not robust, which I assume alludes to the fact that the GTF files don't have the rRNA gene entries necessary, which is the whole point of this question. When I simply looked at the rRNA entries that are in the gtf for hg38p12, there were very few reads reported, but against my 4 rRNA seq ref, there were 5 thousand times more.

And BTW, since I'm using galaxy, do you know of an rnaseqc wrapper that includes that -BWArRNA flag? I can't find one that includes it.

As far as I understand it uses the genes with the "rRNA" biotype in the Ensembl gene annotation.

The doc says it uses the supplied GTF. Is the documentation wrong or am I misinterpreting it?

I just double-checked, and while there are a number of 5.8S entries in the ensemble rRNA biotype annotations for human, none are for the 5S, 18S, or 28S subunits, which is the problem I cited in my question. I tried out SILVA to try and find those rRNAs, not only for human, but for all species.

You appear to misunderstand the purpose of SILVA. It is "a comprehensive on-line resource for quality checked and aligned ribosomal RNA sequence data." It contains multiple rRNA sequences that have been assigned to the human taxon. As far as I am aware, it doesn't contain consensus sequences for any taxa.

And that’s exactly what my question was - whether I was misunderstanding SILVA and if I could use it to count rRNA sequences in an RNA-Seq run. The answer is yes, I was misunderstanding it and no, you can’t use it to precisely quantify specific species’ rRNA contamination.

The SILVA NR99 data only contains 16S/18S sequences. If you want 23S/28S sequences, then try the LSU data from release 132:
https://www.arb-silva.de/no_cache/download/archive/release_132/Exports/

SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits.  While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent.  Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.
You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol_rrna[PROP]".  This returns 240,240 coding sequences (or 30,775 gene sequences).  You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.
The following e-Utilities commands from ncbi will download the gene sequences:
esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta

This is not a common analysis.  Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination.  While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies.  So you have to be careful to count unique reads and not mappings.

SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S.  You can download them from GitHub

I did see that there is a galaxy tool wrapper for SortMeRNA that can filter rRNAs. I haven’t tried it yet. I also had quickly looked through rfam, but I didn’t spend time looking to see what it had for download. I assume the tool doesn’t include the data? I’d be interested to see how the SILVA/rfam data combo compares to my 30k rDNA download from NCBI nuccore...

There does not appear to be a convenient way of downloading only rRNA entries from RFAM.  They appear to have gene sequences in fasta format for roughly 240k rRNAs, but neither their ftp site nor browsing interface provide a way of compiling and downloading them specifically.

I'd recommend using kraken2 for this. There's a pre-built kraken2 database for SILVA (16S only) which can be downloaded from here:
https://ccb.jhu.edu/software/kraken2/downloads.shtml
For other ribosomal subunits, it's necessary to create a kraken2 database from scratch. 16S/18S subunits are in the SSU database, while 23S/28S are in the LSU database:
https://www.arb-silva.de/download/arb-files/
Note that SILVA is a database of collected rRNA sequences. As far as I'm aware there's no canonical "human" 18S or 28S sequence in SILVA, but it includes a collection of multiple 18S & 28S sequences that have been annotated with the human taxon.
Kraken2 has a command-line interface for creating reports:
kraken2 --db <location_of_database> --gzip-compressed input.fastq.gz --report output.kraken.report.txt > output.kraken.txt

Pavian can be used for visualising the results and creating proportion tables.

But SILVA doesn’t have the 5S, 18S, or 28S subunits in it. Any tool suggestion that uses the SILVA database is going to have the same problem. Besides, the question is about what databases of rRNAs to use, not about tools. I can use any tool to count rRNA reads and do other QC, but regardless of the tool: garbage in, garbage out. A tool can’t fill in missing data and SILVA simply doesn’t have all the human rDNA sequences.

Does what? I’ve downloaded all the SULVA Fastq dukes and blasted all 4 subunits against it. They’re not there. I don’t think they’re intended to be there. If they’re there, give me the accessions.

Autocorrect. “Silva fasta files”. Not SULVA fastq dukes”.

No, I mean, show me the SILVA accession number for the specific human 5S, 18S, or 28S ribosomal subunits. They’re not there. I have the sequences for them. I blasted the SILVA fasta files with them. They’re simply not there. And when I use the SILVA fasta file as a reference, I do not get as many reads mapping to that reference as I do when I use the 5S, 5.8S, 18S, and 28S human ribosomal subunits that I manually retrieved from GenBank as the reference.


text:
Polygenic Risk Scores
I have run a GWAS and used PRSice (https://www.prsice.info/) to calculate polygenic risk scores on the data collected from the GWAS. I need to check whether alcohol dependence can predict gambling disorder, through the PRS. However I am having trouble with how exactly to do that. I have a GWAS on alcohol disorder from the one population and a GWAS on gambling and alcohol disorder from the other population, then I calculated PRS for the three GWAS datasets, but how do I use the one to predict in the other? Through the literature it was suggested to use linear mixed models, which I have done but I still feel as if I am doing it wrong as I only did the model on the one dataset and not the two combined.

Before doing anything else, you should split your dataset into training and test groups, and [this is important] don't touch the test dataset until you're fully confident that your approach is correct: no post-test modification.
A split that's considered acceptable in Computer Science is a 10%/90% split, i.e. use 10% of the data for training / algorithm generation and the remaining 90% for testing. This will give you good confidence that the algorithm works.
I personally think it's a good idea to do population sub-sampling (i.e. further subsample that 10%) to determine an informative marker set prior to any risk calculations. What you should be looking for are markers that are consistently linked with the trait of interest (i.e. in your case alcohol dependence) in multiple sub-sampled populations. I go into a lot of detail about how to do this here [bearing in mind it's an old pre-print, and doesn't do the 10/90 split I've recommended above].
After you're happy with the marker set, calculate the risk scores only using that set (hopefully PRSice can do this), then carry that risk score [only] through to test against the dependent trait of interest in the training dataset (i.e. in your case gambling disorder). If your dependent trait is a binary trait, this could be done using ROC analysis (as in my preprint) to calculate area under the curve. Or use linear mixed models, whatever works for you.
Make any changes you think you need to improve the test result (e.g. changing the composition of the marker set).
After convincing yourself that the approach is correct, do the same risk score calculation (using the same set of markers) on the test set, find out how well it predicts the dependent trait, and report on the results. Don't change anything.


text:
When running ChIPQC bioconductor package after peak calling, should I use sorted BAM files or unsorted BAM files?
I have called peaks using MACS2 and now I'm trying to run ChIPQC bioconductor package. Can someone please tell me should I use sorted BAM files or unsorted BAM files? I have created sorted BAM files and from that indexes for BAM files. When I was trying run ChIPQC earlier, I got an error. Then I figured out that is because I don't have index files. Then I created indexes after sorting BAM files.

You need to provide sorted bam files, ideally with an index. If you run ChIPQCsample, it calls ChIPQC:::sampleQC which iterates through the chromosomes or contigs, using the index, and the absence of the index will crash it.
We can use an example, first in linux:
wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeHaibTfbs/wgEncodeHaibTfbsA549RxlchV0416101Etoh02AlnRep1.bam -O Etoh02AlnRep1.bam

samtools index ./Etoh02AlnRep1.bam
# s
samtools view -uh Etoh02AlnRep1.bam chr1 > sorted.bam
samtools index sorted.bam
# quick way to unsort
samtools sort -n -o unsorted.bam sorted.bam
# it fails to index
samtools index unsorted.bam
[E::hts_idx_push] Unsorted positions on sequence #1: 246814554 followed by 42922635
samtools index: failed to create index for "unsorted.bam"

Then in R:
library(ChIPQC)
ChIPQCsample("unsorted.bam")
Bam file has 25 contigs
Error in value[[3L]](cond) :
  failed to open BamFile: failed to load BAM index
  file: unsorted.bam

ChIPQCsample("sorted.bam",annotation=NULL)
Bam file has 25 contigs
Calculating coverage histogram for chr1

Calculating SSD for chr1

Calculating unique positions per strand for chr1

Calculating shift for chr1

 300 / 300
Counting reads in features for chr1

....


text:
DiffBind, diffferentially binding site
I have data for 3 histone marks (2 for silencing and 1 for activation) each mark has three replicates.
when I run the diffBind package I have three contrast:
Silence1-Activation
Silence1-Silence2
Silence2-Activation

When I run dba.plotBox for any of contrasts I will have the same distribution of reads.like below:

I even tried to start with 1000 peaks of each one that showed higher fold change but again got the same pattern for plotbox. How can I see which peaks belong to which histone marks?

most likely the data frame is wrangled wrongly. Even the outliers are the same. This is almost impossible

I got the answer from Rory Stark. the samples should have a unique ID. (https://support.bioconductor.org/p/132502/#132514)

Thank you for your feedback @Mariam


text:
Separating peaks of chip-seq with specific length
My data contain several chip-seq results. I have the peaks called by MACS2.I wanted to only look at those peaks that their size is e.g 500bp to 1000bp. How can I separate those peaks efficiently?
I know one way is to take the peaks.xls file sort them based on the length column and take out the peaks less than 1kb. But I was looking for a more efficient way to do that. Any suggestion?
Many Thanks

I would just write an R script or a Python to read the peaks and write those lines to new files according to their lengths.

You can do it by awk.
awk -v OFS="\t" '{ if ( ($3-$2) >= 500 && ($3-$2) <= 1000) { print } }' in.bed > out.bed

Assuming there is no header.

You can use rtracklayer from R bioconductor to import it, filter and do much more stuff, i use the code from this blog for importing the narrow bed file. Below I use an example file since you did not provide yours:
library(rtracklayer)
fl = "zas1_default_peaks.narrowPeak"

extraCols_narrowPeak <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(fl, format = "BED",
                        extraCols = extraCols_narrowPeak)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |                name     score signalValue
         <Rle>     <IRanges>  <Rle> |         <character> <numeric>   <numeric>
  [1] AB325691     8982-9219      * | zas1_default_peak_1        18     1.31697
  [2] AB325691   16685-18050      * | zas1_default_peak_2        59     1.51283
  [3]        I 101971-102188      * | zas1_default_peak_3        17     1.23717
  [4]        I 102390-103538      * | zas1_default_peak_4        19     1.24967
  [5]        I 152286-152496      * | zas1_default_peak_5        23     1.37631
  [6]        I 175025-175577      * | zas1_default_peak_6        69     1.62063
         pValue    qValue      peak
      <numeric> <numeric> <integer>
  [1]   3.37358   1.89101       146
  [2]   8.12841   5.92232       419
  [3]   3.26012   1.79545        64
  [4]   3.47875   1.97673       708
  [5]   3.89921   2.32948       136
  [6]   9.34429    6.9608       374

The narrowBed and peaks.xls contains the same information. If you really really need to read in the xls, you can do:
library(GenomicRanges)
peaks = read.table("zas1_default_peaks.xls",comment="#",header=TRUE)
gr = makeGRangesFromDataFrame(peaks,keep.extra.columns=TRUE)

Once you have the genomic ranges object, filtering is the same:
 gr[width(gr)>500 & width(gr)<1000]
GRanges object with 208 ranges and 7 metadata columns:
        seqnames          ranges strand |    length abs_summit    pileup
           <Rle>       <IRanges>  <Rle> | <integer>  <integer> <numeric>
    [1]        I   175025-175577      * |       553     175399    186.99
    [2]        I   186807-187489      * |       683     187009    208.66
    [3]        I   339070-339700      * |       631     339258    222.91
    [4]        I   431557-432088      * |       532     431734    195.66

Thank you so much for the answer. it works fine with narrow peaks and bed files but does not work with xls files. import function supports other formats like GFF,TwoBit, Wig, bedGRaph but not xls.

it contains the same information. I have edited my answer to show how u can do the import

Many Thanks, Is there any way to export the files as .xls to further analyze them using diffBind?


text:
data visualization RNAseq : scaling data for PCA and cluster dendogram
I have count data from a RNAseq experiment (2 samples are from normal cells and 3 samples are cells with a disease), and the data is already standardized by trimmed mean of M values (TMM).
I want to do some plots: biplot of Principal Component Analysis (PCA) and a cluster dendrogram to see if the samples normal/disease are well separated (there is a clear difference between them).
Since the data is already standardized (TMM) should I scale and center the data prior to perform PCA and cluster dendrogram??
thank you!

by TMM you meaning running limma or edgeR on it i supposed? Or obtaining the cpms?

Using kallisto I get the TMM, and before looking for differentially expressed genes with edgeR I am doing PCA and cluster dendogram to look into the data.

Scaling (or centering) makes the genes comparable: Putting the expression levels of genes in the same scale (i.e. between 0 and 1) sustains that all of your genes contribute equally to the PCA or distance calculations. On the other hand, without this step, such calculations would be dominated by the highly/lowly expressed genes. I get "better" results with scaling.

Yes, I understand that. But the data has already been standardized using by trimmed mean of the value. And when performing PCA with and without scaling the results are pretty similar. I am wondering if the scaling and centering are not needed because the data has already been standardized.

What TMM does is to "make different samples comparable" by adjusting for library size and more. It does not deal with individual genes. So yes, scaling would help when expression levels of different genes vary a lot (it is often the case). Since you are using edgeR, you can use the function cpm() to calculate log2 of CPM and then scale these before clustering for example.

ok, thank you for your answer now I undestand. But, in terms of using the cpm() function: is it ok to use it over the data that is already TMM? Thanks


text:
How does FDRtool work?
I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = T)
ddsRes[,"padj"]  <- p.adjust(FDR.ddsRes$pval, method = "BH")

Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?
Thanks!!

Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.
Basically, the p.adjust(..method = "BH") is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.
When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.
However, I have to note that in your code, you are not using fdrtool to estimate the fdr as explained above, because:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = FALSE)

With this, the raw p-values are stored in FDR.ddsRes$pval and the q value, which is the corrected p-value, is stored in FDR.ddsRes$qval. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:
p.adjust(FDR.ddsRes$pval, method = "BH")

You should apply a two tailed test. From the wald stats, it should be:
derived_p = 2*pnorm(-abs(ddsRes$stat))

In your code, you are running it wrongly by passing the test statistic into fdrtools, we can compare the raw p-values:
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.


text:
Error in if (is.na(sra_acc$run[i])) { : argument is of length zero
I am trying to write a code for importing sra files as fastq files in Rstudio (for obtaining and aligning RNA-seq reads from public repositories). I keep getting "Error in if (is.na(sra_acc$run[i])) { : argument is of length zero" right after I run this code. Any suggestions on what might be going wrong? My code is as follows:
library(SRAdb)
sqlfile <-'SRAmetadb.sqlite'
if(!file.exists('SRAmetadb.sqlite')) sqlfile <<- getSRAdbFile()
sra_con <- dbConnect(SQLite(),sqlfile)
sraInf <- getSRAinfo("SRP045534",sra_con, sraType="sra")

Error in if (is.na(sra_acc$run[i])) { : argument is of length zero

Have you confirmed whether file.exists('SRAmetadb.sqlite') is TRUE or FALSE? Or generally, what the value of any of your variables is? I am not familiar with this library but having sqlfile be a string that is overwritten by the return value of getSRAdbFile() (don't know what that is) seems like it might cause an issue. Also, are you intending to use the superassignment operator (<<-)? I believe that you are using it in a non-standard way: https://stat.ethz.ch/pipermail/r-help/2011-April/275905.html

After checking with file.exists(), I would suggest you try the same command with another SRA run. It may happen that some runs are not accessible for reasons I do not know. I always have problems in accessing sra formatted files and use a workaround to download fastq.
Here is my not official workaround:
runlist<-sraConvert("SRP045534",sra_con=sra_con)
runlist[1:2,]

study submission    sample experiment        run
1 SRP045534  SRA179382 SRS685546  SRX681992 SRR1552451
2 SRP045534  SRA179382 SRS685544  SRX681990 SRR1552449

for (dloadrun in runlist$run)
{
    getSRAfile(dloadrun,sra_con=sra_con,fileType="fastq")
}

This should iteratively download all the runs associated to the project you are studying.
It is possible that you can do the same with sra type file, but I have never been able to use them with sraDB, and never investigated the issue.

I am also having this issue- my SQL file comes back as TRUE when I use file.exists but then in Fabio's workaround, the table comes out blank. Same error as you if I use getSRAfile.


text:
How do I create a filtered gene list using expression medians
Forgive the simple noob question
I have TPM data of ~50k genes (rows) across ~1k cell lines (columns). In R, I would like to output an "intermediate expression" gene list for each cell line, like:
>head(finalExpressedGeneList, n=3)

Celltype1    Celltype2    CelltypeN
gene1        gene5        gene3
gene3        gene6        gene6
etc          etc          etc

I am defining intermediate expression as "gene expression >= median value within each cell line." I have already replaced all 0 values with NA using:
>data[data == 0] <- NA

I then got medians for each column using:
>colMed <- apply(data,2, FUN = median, na.rm = TRUE)

Now I need to somehow filter the rows from each column and copy the associated gene into a new matrix. I started building the code below:
>if (data[,1]>=colMed [1])
   {
   finalExpressedGeneList <- row.names(data)
   }

But I realized that even if I got it working, I would still have to iterate over all of my columns. I'm sure that there is a cleaner way to do this.
EDIT (adding example data)
#create example data
>j = 5
>M <- matrix(NA,j,j)
>M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
>rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
>colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
>view(M)

       cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
gene1          1           1           1           1           1
gene2          2           2           2           2          NA
gene3          3           3           3          NA          NA
gene4          4           4          NA          NA          NA
gene5          5          NA          NA          NA          NA

I would like to get an output of:

cellLine1   cellLine2   cellLine3   cellLine4   cellLine5
    gene3       gene3       gene2       gene2       gene1
    gene4       gene4       gene3
    gene5

because the median of cellLine1 has three genes greater than or equal to the median of that column (3), cellLine2 has two genes that are greater than or equal to the median of that column (2.5), etc...

This below should work. It returns a list,
## example data
j = 5
M <- matrix(NA,j,j)
M <- sapply(1:j, function(i) `length<-`(1:(j-i+1), j))
rownames(M) <- c("gene1", "gene2", "gene3", "gene4", "gene5")
colnames(M) <- c("cellLine1", "cellLine2", "cellLine3", "cellLine4", "cellLine5")
#####
genelist = apply(M,2,function(i)rownames(M)[which(i>=median(i,na.rm=T))])

####
#you can call the genes for cellline 1 using:
genelist$cellLine1

Hope I got it correct!


text:
Determining the quality of a publicly available genome assembly
I am looking for homologous genes in various protozoan species using BLAST. The genome sequences of these species are deposited in NCBI's WGS database. NCBI’s webpage includes the global statistics for each WGS project and it indicates the assembly level (scaffold, contig, etc.), the assembly statistics (number of contigs, N50, L50, etc) and (if available) the genome coverage. Can any / some of these parameters be used to evaluate the genome assembly’s "quality", especially if there is more than one genome assembly for the same species? Can these data be used to determine if the genome is very fragmented or not?
I understand that quality is a relative term and that there are several options for assessing genome assemblies (QUAST, BUSCO, etc.) but I am just looking for a "simple" way to compare several genome assemblies (for e.g., this post How can I improve a long-read assembly with a repetitive genome? states that the reference genome is highly fragmented, but there is no explanation regarding how the OP determined this). I believe that the N50 value can be helpful, but I need some context to understand what is "good" or "bad".  I am at a loss of where to start, so any information will be greatly appreciated!
-Leah

It is not so easy to determine the quality of an assembly - but sequence length, number of contigs, N50 and L50 can give you some guidance.
How much does the sequence length deviate from the expected length (or how much variation is among assemblies of the same species). Same for number of contigs/scaffolds. The sequence length, and number of contigs together with N50 and L50 give you an idea about the fragmentation. Very small N50 compared to expected genome size and chromosomes is a likely indication of high fragmentation, same for high L50 compared to expected number of chromosomes.
However, these number can be misleading. High N50 and low L50 can occur witho some very long contigs, so you may want to compare N75/L75 or better the Quast statistics

Carambakaracho's answer is good.  A high N50 should usually be a good sign, as well as agreement between the overall genome assembly length and the expected genome size.  The number of 'N's in the assembly may also be important to check, as Ns are sometimes used to indicate single unknown bases, and sometimes to indicate unknown numbers of unknown bases, and Ns may be used to connect scaffolds,which may lead to artificially higher N50.  If you have access to a set of genome sequence reads (Illumina) for the same organism, then another more detailed check you can make is to compare the kmer spectrum of the reads with the kmer spectrum of the assembly, using the kat tool.
Finally, if you want to compare several assemblies, you might align them against each other, to see how they compare qualitatively (e.g. using mauve).

Thank you very much for your answers, I really appreciate your help. @Jonathan Moore, just a couple of questions: (1) what would be considered a good / bad kmer spectrum of the reads? and (2) is it possible to align the whole xxx_genomic.fna files from NCBI using the mauve tool?

a good assembly should have similar kmer spectrum for the reads as for the assembly, with the ratio in abundances between the two spectra reflecting the depth of coverage of the reads.  If some kmers have a large difference in abundances in the reads compared to the assembly (double for instance), it may be that these parts of the assembly do not have haploid copies collapsed correctly, or some other issue with the assembly, and if there are many kmers in the reads which are present in large numbers but absent from the assembly, it may be evidence of missing parts of the genome.

I think mauve should be able to align the whole genome .fna files. It may be somewhat slow with larger genomes.

Thank you very much @Jonathan Moore! I will give them a try.


text:
How to make a robust amino acid phylogeny with adequate robustness support
This is the second part of this thread:  Duplicated genes and genome fragmentation
As @Michael G. suggested, I did a pseudo-annotation of sequences 5-2 and 5-3 (i.e., I joined both sequences into one) and I did the phylogenetic analysis again. Here is the tree that I obtained:

Both of my sequences are shown in a red rectangle. The length of the branch seems to have decreased a bit. Also, I looked at the branches near sequence 5-1. One of them is the original query sequence from species A and the other one is a hypothetical protein from the same species (probably related to the query sequence, but I couldn't find any domains because it is only a few amino acids long). I still need to find the root of the phylogeny but I am not sure how to do this...

So what is the actual question here? How you can root the tree?

The question is how to make a robust amino acid phylogeny with adequate robustness support. This is a generic question that many people will want answering and something I know a lot about. In my opinion, this has been a successful investigation given that the final tree is not a "rough and ready tree", but a formal likelihood-based tree. I see this as the end of the investigation and final response

I just have not got around to answering at this moment. I hope this is ok and we'll know in future when a rolling investigation is in-play either we  section the questions to specific topics or find an alternative forum

So @leah for forum rules it might be better to delete this and ask again OR change the title and specifically what you need is state of the art tree building to finish and conclude your study. Amino phylogeny is complicated and not trivial so in my opinion is a valid question. The idea is that someone in years in the future can understand from this post how to address a specific question. Stack exchange aims to get 90% of traffic via Google questions. Thus the question format is off-topic.

Exactly, thanks @MichaelG.! I agree that you two have had a very interesting discussion, it's really a shame it happened here since this site is not a good place to store discussions. It would be really, really great if you could combine all of this into one comprehensive Q&A though since there is a lot of great info you've both come up with.

Ok.. @leah I will try and set up a Q&A, I can't remember how to do this, but I'll give it a whirl. Essentially the site has a chat feature.

Thank you for pointing out the problem with this thread. As @Michael G. suggested, I changed the title so that it is better suited to the format of the forum. Hopefully these posts can be moved to the Q&A format so that the information is available for other people later on.  As for the question itself, I agree with Michael G. that we are almost done, and I am very much looking forward to his input to round things off and get some final conclusions on this investigation!  Again, many thanks for all your help, I have certainly learned a lot through these discussions.

I cannot comment... which seems odd while I can write answers.
Could you share with me the data you are using for the tree? I am developing a software for phylogenetic tree generation. I'd like to compare the tree and see if it can help you.
You can also try yourself if you like: https://www.dnaservic.es/
I find that the root it finds makes usually far more sense than with other algorithms.

Thank you @juanjo75es. I will give your software a try and I will let you know how it goes. I hope the results will helpwith the software you are working on.

To finish you need bootstrap support using a model-based tree building algorhithm, I use RAxML, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.
A robust maximum likelihood tree is,
raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,
raxmlHPC -b 459 -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

Make a consensus tree of the bootstraps,
raxmlHPC -f b -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)
You need a fairly robust desktop, dedicated for around 24 hours to complete one calculation. You have enough results to draw a conclusion.
Viewing the tree, FigTree is easy.
Root
Check the work of Cavalier-Smith, he heavily focused on basal protists to find a species which evolved, using external evidence, before any other member of your group. One of the ingroups may have sufficient evidence to say it evolved before the rest and that is your root, which you select as an outgroup.

The mod is asking for a Q&A, when I find out what this is, this will need to be done.

Thank your for the detailed instructions @Michael G. I will give it a try and let you know how it goes. And thank you for the heads up about the program taking a long time to complete the calculation... I was not expecting that!

Hi @Michael G. Just a quick note to let you know that I am still working on the phylogenetic tree, but since I have no experience compiling, it is taking me a little longer than I expected. I will post the results as soon as possible.

Hi @leah you can always post it as a separate question

That's a good idea @Michael G. Will do!

Well, it has taken me some time but I have finally managed to install RAxML and FigTree in Linux (I could not compile RAxML on OSX). I took my alignment file, saved it in Phylip format and I tried to run the first command line to build the robust maximum likelihood tree (raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s ~/Desktop/alignment_file.phy), but I got the following error: Error, the outgroup name "outgroup-p" you specified can not be found in the alignment, exiting ....  and I don't know how to fix it. Could you please help me out? Thank you very much!!

Can you put this as a separate question please?

Certainly. The question is posted in this thread: https://bioinformatics.stackexchange.com/questions/10534/how-to-define-an-outgroup-to-build-a-robust-amino-acid-tree. Thank you @Michael G.