text:
Is it conflicting: significant difference according to Wald test but not in PERMANOVA and ANOSIM to compare species abundances of metagenomic samples?
If there is no significant difference between two groups of metagenomic samples according to "Multivariate differential abundance tests" for example PERMANOVA and ANOSIM, does it makes sense to use "Univariate differential abundance testing" to see, for example, what species have significantly different abundances between groups?
The issue is that I found species that have significantly different abundances among samples of two groups using Wald test in DESeq2, but there is no significant difference between the two groups according to PERMANOVA and ANOSIM.

Yes, it does make sense to perform differential abundance testing after one doesn't find compositional differences with PERMANOVA / ANOSIM.
It seems to me you are drawing an (implicit) analogy to ANOVA and post-hoc tests like the Tukey Honest Significant Difference test. Intuitively, the analogy makes sense, however, these situations are quite different. PERMANOVA / ANOSIM and DESeq2 are operating on different aspects of the data, with different transformations and under different assumptions.
PERMANOVA and ANOSIM do not operate on the raw counts directly, rather, they test differences between groups by calculating some distance or dissimilarity between the samples. One can choose between several distance / dissimilarity measures, and the test result may be impacted by the chosen measure.
DESeq2 (and other packages), on the other hand, test for species abundances under a parametric model, operating on the raw counts directly.
Thus (sadly) one is not a natural, optional follow-up of the other, they are complementary and will tell different stories about the data.


text:
parsimony and maximum likelihood tree comparison in R
Maximum Likelihood analysis compares the tips or the species based on their sequence similarity while Parsimony analysis compares the characteristic features among the species.
So far I did comparison of both the tree using dendextend. Is there any other methods or library to make the comparison more quantitative.Due to my limited understanding I have this in mind. Can i even compare these two methods?

Yes you can definitely compare both methods by tree topology, personally I do this via boostrapping.

as usual what is the R way of doing it?

How many taxa are there in your tree?

this is my parismony data "Phylogenetic tree with 50 tips and 48 internal nodes." and this is my ML  Phylogenetic tree with 50 tips and 49 internal nodes.

@Michael shall i share both the tree? object here?

Hi no that is fine ... it is just that ML takes a long time and the number of taxa is key to the computational effort .. okay I will look at phangorn and post the code to bootstrap it. Do you use R on a cluster and are you happy 'farming' a job, e.g. splitting it in two or three?

"Do you use R on a cluster and are you happy 'farming' a job, e.g. splitting it in two or three?" no i use workstation but we do have a cluster..i can use ..

Ok I'll post up the codes, I need to look at phangorn it certainly does bootstrapping, but mid-week if thats ok? You know deadlines, deadlines ... :-(

I don't use R to build trees, I use established subprocessing stripts in Python and ETE3. So I've not run this code, but I do understand the critical parts of the code.
The hard bit is bootstrapping a maximum likelihood tree, first is the ML model
1.
fitGTR <-optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE, + rearrangement = "NNI", control = pml.control(trace = 0))

The above model is perfect, general time reversible model with a gamma correction and an estimate of the invariant sites. This is the most important part of tree estimation. The NNI for tree rearrangements is a good one. So THIS is the model you want. You can emperical test this, but these days we don't bother.
The information given for bootstrapping,
2.
BStree <- bootstrap.pml(x, bs = 500, trees = TRUE, multicore = FALSE,
  mc.cores = NULL, ...)

The critical parameter is bs=500. By default it uses bs=100 which I don't think is enough. I would set multicore = TRUE to speed the calculation (it will speed it up alot) and I assume there must be a setting for the number of cores (?).
I assume you can then stack all the parameters from 1/ into code 2. So in other words the model you used in part 1. (GTR PINVAR and Gammma) is same model used in the bootstrapping.
This will produce a bootstrap file, probably of 500 bootstrapped ML trees, so you need to produce a consensus and represent the values on the optimal ML tree.
Thus you need to overlay the bootstrap onto the tree that is output in 1. This says fitGTR and I suspect that will be a tree file with a likelihood, you might need to set a tree output and that would be 'tree' in the code below.
3.
plotBS(tree, BStrees, type = "myroot", bs.col = "black",
  bs.adj = NULL, p = 50, frame = "none", ...)

The type = myroot is your defined outgroup for the tree. This should be your nicely bootstrapped tree.
There was also this code in the help manual and I'm not sure what this is doing. I suspect it is just producing bootstrapped datasets, so its not needed here.
bootstrap.phyDat(x, FUN, bs = 100, multicore = FALSE,
      mc.cores = NULL, jumble = TRUE, ...)

To do the same with parsimony then just replace pml for parsimony. Comparing 2 bootstrap trees is a node is robust when it scores >80%

will use this one to run

Can i even compare these two methods?

See if this paper helps:
https://www.nature.com/articles/nature02917

Welcome  @HD2000 thanks for your comment, could you please summarise the content of your weblink here? I agrree it is a very thought provoking paper indeed, one that didn't stand the test of conventional practice but waooo

Can you add more detail about why your link would help?

source: Stack Exchange

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

sort: votes


text:
What unit (TPM, FPKM/RPKM, or other) to use when working across samples
I have raw gene read counts and would like to perform an analysis across multiple samples. I've found conflicting info online on how this should be done. One commonality however is that FPKM/RPKM should not be used. On the other hand, there seems (to me at least) to be conflict about whether TPM is suitable. For example, this StatQuest video seems to suggest that TPM is suitable for analysis across samples. One person on this Biostar post (among other resources I've looked at) suggested that none of the units (FPKM, RPKM, TPM) are suitable for cross sample comparisons.
I am quite confused. Can someone clear up the confusion?

Choices in bioinformatics tend to be quite situation dependent; could you please provide more information about the specific problem you're having? What tool are you using for mapping? Are the reads long reads, paired end, or single reads?

Neither of RPKM, FPKM or TPM is a good choice. All these methods are similar, they only perform the operations in a slightly different order. Eventually they correct for gene length and sequencing depth. Correction for sequencing depth is necessary but not sufficient in many cases. Correction for length is typically not done for a differential analysis as it reduces the counts of many genes and therefore throws away statistical power. I will not describe this in detail since the StatQuest video series you link includes videos about DESeq2 and edgeR normalization procedures which extensively cover the normalization procedures and why these approaches are superior to RPKM/FPKM/TPM. In short: You have to correct for differences in library composition to compare between samples. Please watch these videos and read the manuals of the established RNA-seq analysis packages such as but not exclusively limma, edgeR and DESeq2.

And what about EBseq?

What about it? You do not mention it in your post, please explain.

I mean could you use it similarly to edgeR and DESeq2?

Never heard of the method. Probably you can, I suggest you read benchmarking papers to learn how it performs and whether people recommend it.

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


text:
FASTQ file trimming
I have a FASTQ file:
@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
ACTTANGTGTATGTAAACTTCCGACTTCAACTGTATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
ACTTANGTGATGTAAACTTCCGACTTCAACTGTATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

I want to remove the string  ACTTAAGTGTATGTAAACTTCCGACTTCAACTG from it while retaining the sequences starting from TA onwards.
I wrote grep "ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA" SRR_1.fastq | sed "s/ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/g", but it removes the headers of the FASTQ file as well.
The desired output would be
@NS500455:80:HG7TNBGXB:1:11101:17723:1055 1:N:0:ATCACG
TATAGGGATCCNAGCTCCAATTCGCCCTATAGTGAGTCGTAT
+
/AAAA#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE#EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@NS500455:80:HG7TNBGXB:1:11101:8821:1057 1:N:0:ATCACG
TATAATAAATATCTAGATCGGAAGAGCACACGTCCGAACTCCA

How to do this. I dont want to use cutadapt.
Kindly help.

What application are you using that will be okay with a mismatching quality string?
Use trimmomatic or cutadapt, or bcl2fastq could probably do the job too.  Those will also have the benefit of dealing properly with one-off errors, which awk and sed won't.

Skip the grep, use just the sed. Unless the header contains the long sequence as well (which I doubt it will), your command will work.
sed 's/^ACTTAAGTGTATGTAAACTTCCGACTTCAACTGTA/TA/' SRR_1.fastq > SRR_1.edited.fastq

It works but there are some sequences begining with: ACTTANGTGTATGTAAACTTCCGACTTCAACTG They are being retained. How to deal with them. Kindly guide.

That's because N does not match A,T,C or G as a literal. Replace all N with [ATCG] to remove these sequences.

It is worth considering thomas' point. If you'd like your QUAL lines to reflect the change to the SEQ lines. awk or bioawk might be better at this point.

If you need the quality vector for each read, then your code will create a shift because it does not cut the quality vector. You should rather use a trimming tool a systematically cut all 33 first characters. fastp will allow you to do it easily with the -f, --trim_front1 argument.

What should be the command line , I typed fastp -i SRR_1.fastq -o trimSRR1.fq -f 33 and it gives an error ERROR: trim_front1 (--trim_front1) should be 0 ~ 30, suggest 0 ~ 4. Kindly help.

strange! I did'nt know it was restricted to 30. You can run fastp twice with 20 in the first trimming and 13 in second one? not really convenient but it should work!

you may alternatively had a look at Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). It has a HEADCROP argument, but I don't know if it also restricted in length.


text:
minimum number of neighbors for scRNA-seq analysis
Many scRNA-seq clustering/visualization methods are based on nearest neighbors. The maximum number of neighbors is limited by the neighborhood size. Is there a method to determine the minimum number of neighbors? Now that most datasets are well in the thousands of cells, this point is largely irrelevant, so maybe that is why not much guidance on this issue is available. However, if you something like 50 cells, it becomes far more critical.
If you have a deeply sequenced dataset, your neighbors should be more reliable (less noisy) than one with shallow coverage, so you should be able to use less. Is there a way to somehow quantify that?

Honestly I feel like you could do a PhD around this question.


text:
Clarification about how mapping to a reference genome works in RNA seq
When mapping reads to a reference genome, how is it possible to tell which part of the genome the read is referencing if the genome possibly has the same repeated sequences in different areas? What I mean by this is let's say for the sake of argument that you have a read "AGCTAGCT". If the reference genome has two places where the sequence "AGCTAGCT" is present, how do we tell which one it's coming from? I'm assuming that the length of the read is chosen long enough such that this doesn't happen but perhaps I'm wrong?

Of course longer reads will decrease the chance of having such a situation. If you use Bowtie2, you have a randomized search for equal mappings (http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml, see the Randomness in Bowtie 2 in the manual). In addition, it will report reads that align, 0, 1 or more times.

@thomasdugedebernonville I'd put that as an answer rather than a comment.

As the other answers/comments have covered, some reads will inevitably align to multiple positions in the genome.
This can come from a variety of sources, like aligning to repetitive regions of the genome, short read fragments (that will likely match multiple regions in the genome, even if it, in theory, came from a single location), and low base call qualities (the aligner doesn't have lots of confidence aligning it it to one location because it could map to another, with some likely mismatched calls in the read).
This isn't inherently a problem, depending on the analysis you want to do after aligning.
In fact, having a single read that is compatible with multiple transcripts is one of the ideas behind the reference index used in kallisto, and the "transcript compatibility counts".
If you want to minimize this effect, two easy ways are to create longer reads on your sequencer, or to use paired-end reads (in theory, both ends should align relatively close to each other, so if one of the mates maps to multiple places, you can borrow from the confidence of the alignment of its mate).

You can't always know exactly where every read originated, even if you use long reads to sequence full-length transcripts (odds are there's at least one transcript that's duplicated in the genome). This is also not exactly a problem, such reads can either be discarded or reported as multimapping.


text:
Gffcompare issue: 0 reference transcripts loaded
I am trying to use gffcompare to compare my assembled transcriptome to a reference gtf that contains information about small open reading frames (sORFs). The reference gtf was obtained by processing downloaded data from several sORF databases and is called all_38.gtf.
I used the command: gffcompare -r path/to/all_38.gtf -o /output_folder my_transcriptome.gtf.
mytranscriptome.gtf is generated by assembling my Bam files with reference to ensembl hg38 v99 reference.
Upon completion, the output of the error message log file reads:
0 reference transcripts loaded.
  237788 query transfrags loaded.
  2714 duplicate query transfrags discarded.

The expected output files are generated except for the .refmap file which I need for downstream analysis. I assume this is because 0 reference transcripts were loaded into the program. Has anyone encountered similar issues? And is possibly something wrong with my custom reference gtf file?
My reference gtf can be found here. The first few lines of my input gtf look like this:
   # StringTie version 2.1.2
1       StringTie       transcript      11869   14409   1000    +       .       gene_id "MSTRG.1"; transcript_id "ENST00000456328"; gene_name "DDX11L1"; ref
_gene_id "ENSG00000223972";
1       StringTie       exon    11869   12227   1000    +       .       gene_id "MSTRG.1"; transcript_id "ENST00000456328"; exon_number "1"; gene_name "DDX1
1L1"; ref_gene_id "ENSG00000223972";
1       StringTie       exon    12613   12721   1000    +       .       gene_id "MSTRG.1"; transcript_id "ENST00000456328"; exon_number "2"; gene_name "DDX1
1L1"; ref_gene_id "ENSG00000223972";
1       StringTie       exon    13221   14409   1000    +       .       gene_id "MSTRG.1"; transcript_id "ENST00000456328"; exon_number "3"; gene_name "DDX1
1L1"; ref_gene_id "ENSG00000223972";

EDIT: Output when running with --debug option
Running with debug option set results in an error file with more information about duplicated transfrags but not about reference transcripts. I've included a few lines of the error file:
Prefix for output files: norfs
Loading reference transcripts..
  0 reference transcripts loaded.
  237788 query transfrags loaded.
        Query transcript ENST00000655384 discarded (duplicate of ENST00000667728)
        Query transcript ENST00000415295 discarded (duplicate of ENST00000670824)
        Query transcript ENST00000263741 discarded (duplicate of ENST00000660930)
        Query transcript ENST00000473215 discarded (duplicate of ENST00000347370)

The --debug option also generates 3 additional files:<outprefix>.Q_discarded.lst, <outprefix>.missed_introns.gff, <outprefix>.R_missed.lst.
However, in my output folder, I only have the missed_introns.gff and R_missed.gff files. The Q_discard.lst is missing.
One thing that I'm not sure matters is that my_transcriptome.gtf was generated by aligning and assembling with reference to an Ensembl gtf file. The gtf that I am using as reference input in gffcompare is downloaded from norfs.org. It contains information about small Open Reading Frames (sORFs). Could the different gtf files used be causing this issue?
If so, the solution would be to assemble my transcripts using the sORFs gtf as reference. However, I'm not sure if this is a good idea since sORFs have relatively short lengths (< 300bp) and this might screw up the assembly.


text:
Regarding RNA seq data analysis and building coexpression network
I have some questions regarding RNA seq analysis if you can suggest anything it will help me a lot.

I am currently normalizing RNA seq data for comparing genes expression within and between samples. Which normalization method would you recommend for this type of analysis? FPKM, TPM, TMM? Also I want to make a heatmap to see genes expressed in different conditions. Do you think transforming normalized data (like log2, z-score) is a good idea for this?
Also I want to build a co expression network, I am just wondering if normalization like FPKM, TPM, TMM  has any influence on building a coexpression network?
Another thing I want to do is use Pearson correlation but I am confused if it will only find linear relationships among the normally distributed data.  But normalization methods do not assume that counts to be normally distributed. So, Is it a bad idea to find pairwise coexpression of genes using pearson correlation? If so, which method would you recommend is reliable for building coexpression networks?

Please help me with this, I will be very grateful to you
Thanks in advance


text:
Which module to select for Pathway analysis based on Module trait correlation and pvalue?
I have a total of 35 tumor samples classified into 4 subtypes. Subtype A, B, C, and D. I have RNAseq data. I'm interested in identifying modules related to each subtype with co-expression network analysis. I have read the WGCNA tutorial and create a trait.csv file like below:
Showing here some rows of the trait.csv file.

SamInfo <- read.csv("trait.csv")

Samples = rownames(datExpr);
traitRows = match(Samples, SamInfo$Samples);
datTraits = SamInfo[traitRows, -1];
rownames(datTraits) = SamInfo[traitRows, 1];

nGenes = ncol(datExpr);
#15648L

nSamples = nrow(datExpr);
#35L

#calculate eigengenes (1st principal component) of modules
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

#Will display correlations and their p-values
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")",
                   sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));

#Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),
               yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE,
               colors = blueWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE,
               cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))

The Correlation between Subtypes and Modules Heatmap looks like below:

For functions identification (Pathway Analysis), from the above output for each subtype, should I select only the modules that were significantly positively correlated with subtypes, or Can I also use the significantly negatively correlated modules?
Is it good to use the modules that were not significant for pathway analysis?
In case if I can use positive and negative correlated modules, Can I combine those modules and use those genes for Pathway analysis?

It's a great question. Please consider the answer with a view to upvoting - or even accepting - the answer. It's good all round.

If you can run pathway analysis programmatically on all modules, yes, do it on all modules. Otherwise, select the modules with strongest associations with each subtype, both positive and negative. I would keep the modules separate because they may have very different biological functions.

Thanks for the answer. One small question. Lets say from the above image, for Subtype A I have selected modules salmon, grey60, tan (negative), and midnightblue (positive) for pathway analysis. I see that salmon module is also significantly positive with Subtype C. So, then the pathways I get from salmon module are related to which subtype?

Well, I don't see the figure anymore so I can only go by your description, but if salmon module is negatively correlated with A, it means (most of) the genes in the module are downregulated in subtype A and upregulated in (positively associated with) subtype C.

thanq. sorry for the heatmap with correlation values modules and trait information. Have some doubt whether here MEs0 = moduleEigengenes(datExpr, ?)$eigengenes I have to use dynamicColors which is before merging or moduleColors which is after merging. In the above image I used dynamicColors. 2) there are few modules like turquoise and blue with more than 8000 and 5000 genes. The number is huge. How do I reduce the number?

You can use either merged or un-merged colors; I would suggest using the merged ones.

Some modules are often very large which just reflects a factor (driving expression or co-expression) that strongly affects a lot of genes. This could be biological or technical, and could be important to you or a nuisance confounder. You could try to identify the factor e.g. from sample clustering (could be a batch effect) or by relating the MEs of largest modules to known sample characteristics (age, sex, race etc). If you identify a match, adjust the data for it, or simply adjust the data for the ME and run WGCNA again.

thanq very much. one more question. To visualize all the modules as network plot....Is it only cytoscape an option? Is there any other software or code to create network plot....


text:
Calculating Synonymous and non-synonymous
library(seqinr)
library(ape)
library(phangorn)

alignment file data file
    sylvia.aln <- read.alignment("abhi_seq/all_seq_no_gal.fas", format="fasta")
    sylvia.aln
kaks(sylvia.aln, verbose=TRUE)

Warning message:
In kaks(sylvia.aln) : sequence lengths are not a multiple of 3

It ends with this warning and there are no result.
When i ran one of the seqinr example data what I see there are no gaps in the sequence files. If i understand the results is due to the gap I suppose.
Any help or suggestion would be really appreciated.

Woaa that is complicated data.
The reason it fails is easy, correcting it is harder. It fails because the open reading frame is lost during in silico translation, e.g. several gaps have 4 indels and the final nucleotides after the final indel are always 4. Thus it will invariably encounter a stop codon OR it will simply refuse to deal with an indel if it is incorporated into a triplet codon. There is no triplet codon model for an indel so the whole calculations stops.
Correcting data that complex ... the easiest way to do it is align you entire data set using amino acids, e.g. muscle. Have an identical nucleotide data set. MAKE SURE the names and positions are identical between your amino acid and nucleotide data, as a programmer this isn't hard to do however, everyone has their own scripts to do this. If you don't do this you next step will fall over.
Use a amino acid to nucleotide alignment translator. I use transalign (I think its EMBOSS). This will give you an inframe alignment. I noted many/all of you start codons are ATG which is a good sign. Check the nucleotide output and make sure there is no overhangs of the triplet codon at the terminal 3' end and check in Jalview or whatever your alignment editor of choice is that the translation is okay.
This will alignment work in the synonymous/non-synonymous analysis and all subsequent/complimentary analyses.
If you have an intron in protein your you'll need to remove it (hopefully not).
Finally you need to be careful of '3rd codon saturation' in your data because that heavily affects your trees and its certainly a risk given the heterogeneity of the indel structure, but first I would get an inframe alignment. There are checks for saturation (I'm sure even R can do this).

"he easiest way to do it is align you entire data set using amino acids, e.g. muscle." this i did in MEGA .Yes I have respective amino acids for the set of nucleotide as well "MAKE SURE the names and positions are identical between your amino acid and nucleotide data"is it the file name?


text:
Is outlier detection the right choice for determining if a value lies significantly away from the mean in a population in gene frequency?
I have the following dataframe that is comprised of enzyme substrates as the variable and then the unique number of times enzymes capiable of degrading this substrates appeared in my genomes and then,total count of times all of those enzymes appeard in my samples. For example, a 7 unique enzymes were found to be capable of degrading chitin, and they occured a total of 298 times across my genomes.
I wish to see if any particular substrate has significantly larger occurence in both unique and total, than the rest of the substrates. Would outlier detection be appropriate for this or would a simple T-test do?

I agree you have the right idea because that data doesn't look at all linear, so regression analysis looks really complicated. Your approach a straight comparison keeps it simple and is within the biological question you are asking,
However, I would use a non-parametric such as Wilcoxon instead of a T-test. I don't fully know the null distribution you are comparing it against. Wilcoxan doesn't require an assumption of normality. T-tests invariably require a transformation to normality and that gets complicated, but a log or sqrt transformation might be okay. If Wilcoxan produces a no significant difference you report it and move onto a T-test.
It is difficult to use outlier analysis with 9(?) observations, thats used if you have a way-out result and can afford to exclude it (>3 standard deviations). Outlier analysis is to remove a single result that is messing up your analysis and doesn't affect the bulk of the data. You alot of observations to be considering outlier analysis.

The actual dataset is 130 substrates thankfully. But I still think you are correct that it would not be normally distributed. Thus would a glm work here with either a negative binomial or poisson as the distribution?

I think negative binomial would be best as I don't think the variance is equal to the mean, but I guess Icould just test that assumption.

That would need fitting and at thiss point you are into formal regression analyses

I have changed the analysis now to have the frequence of occurance of each enzyme and lables them via the enzyme substrate so my model is glm.nb(substrate ~ frequency). I think this is better as it generates a mean and variance for each substrate thus allowing a much better identification of if the degridation of any singular substrate is enriched?

Okay thats a different question


text:
Files for paired end RNA sequencing
I am looking at the videos at a DIY Transcriptomics course and the speaker mentions that to run Kallisto for read alignment with paired end sequencing, one would enter the following:

where the sample files at the end refer to the two different end reads (forward and reverse read files). However, browsing the NCBI for paired end data, I noticed that only 1 Fastq file is provided per replicate for a paired end procedure (here's an example). So I'm a bit confused about this.

Building off of the answer by @An Ignorant Wanderer:
Essentially there are indeed the paired end data files available, but it takes a little bit of extra navigation to reach it.
Using the specific example you provided of BioProject PRJNA574273:
First select a specific sequencing "Run," aka the sequencing of a specific sample / experimental condition. For example, click on the link to run: SRR10189234
This brings you to the default tab: "Metadata." As @An Ignorant Wanderer pointed out, then you can navigate to the "Data Access" tab.
Then look under the "Original format" subheading. Here you'll find links to two files, corresponding to the two orientations of the paired end sequence files.
Edit: Also please see the answer by @burger, the website linked in that answer is much more user friendly in my opinion.

The problem is that the SRA website can be a little confusing. View the same project on ENA where the organization is more human-readable: https://www.ebi.ac.uk/ena/browser/view/PRJNA574273

Actually after clicking on one of the "Runs" and going to "data access" I noticed two files are provided per run

Paired ends is a configuration in sequencing platforms. Illumina is the most popular at the moment and most single-cell sequencing is done using paired-end Illumina sequencing. It varies between single-cell technologies but most 3’ scRNA-Seq protocols use paired ends.
Kallisto is a tool for handling raw data. Many databases supported processed data in BAM format or UMI/read counts. You should not need to run kallisto unless you are involved in handling raw data from sequencing performed by your laboratory. Some databases such as SRA and EGA support raw data in FASTQ format, others do not for ethical reasons or the data is too large.
I recommend downloading from SRA via their command line interface (this is documented on their website). It is possible to extract paired ends from an aligned BAM file using bedtools and samtools. I’ve done this with the Macosko et al (Cell 2015) DropSeq data for example.
Note that in some cases, the barcode is stored in databases as the headers of the FASTQ, BAM, BED, or CTSS file rather than a separate read. Single cell transcriptomics is an emerging field so there isn’t a standard way to do this, it varies between labs and data centres. For single-cell specifically, some databases are being set up, for example: https://data.humancellatlas.org/


text:
merge rows (to columns) in R dataframe based on IDs
I have a dataframe with two columns: $1 = transcript ID, $2 = expression counts (tpm). This comes from a merging of several kallisto counts (if it matters).
I would like to build, in the most basic R way (no tidyverse), a dataframe with as column $1 the unique transcript IDs, and the expression countsas following columns (as many as needed).
Sample input:
$ head test1
G01000001_1 0
G01000001_1 0.00855717
G01000001_1 0.0113467
G01000001_1 0.0121086
G01000001_1 0.0121108
G01000001_1 0.0148402
G01000001_1 0.0183057
G01000001_1 0.0268394
G01000001_1 0.0298547

Sample output:
$ head -3 test2

G01000001_1 0  0.00855717 0.0113467 0.0121086 0.0121108 0.0148402 0.0183057 0.0268394 0.0298547 0.0444686 0.046795 0.0547494 0.0640871 0.0674897 0.0719834 0.082262 0.0867695 0.094905 0.10269 0.11488 0.121549 0.131143 0.141432 0.145758 0.167886 0.184162 0.192757 0.200921 0.249595 0.283296 0.430173 0.473344 0.518981 0.674573 0.701431 0.845816 2.77479 6.28394
G01000002_1 0  0.0277237 0.0294273 0.029438 0.0379258 0.0606153 0.0614727 0.069608 0.142541 0.197644 0.218152 0.22193 0.227401 0.230179 0.259295 0.279463 0.289982 0.292763 0.303268 0.305259 0.308236 0.319988 0.532091 0.632095 0.707766 0.746436 0.821887 0.822496 0.997334 1.06616 1.21752 1.48436 10.0862 2.01763 3.94584 5.02089 7.55767 8.3501
G01000003_1 0  0.00265959 0.00322424 0.00920435 0.454749 0.692183 0.709159 0.948286 1.11827 1.39854 1.39985 1.5781 1.86343 10.2554 10.7185 11.0315 14.342 15.2383 2.13068 2.41792 2.52833 2.77868 3.56981 3.84991 3.97521 4.10445 4.16344 4.26902 4.97202 47.8293 47.9114 5.21789 5.21926 5.35712 5.49942 5.7085 5.87654 50.7484 51.3721 56.7279 7.35062 7.54485 7.64363e-05 8.1367 8.34828

This could be achieved in awk, but I would really like an R solution.
For the logic, here is the awk way:
sort test1 | awk '{if(a!=$1) {a=$1; printf "\n%s%s",$0,FS} else {a=$1;$1="";printf $0 }} END {printf "\n" }' > test2

For a data.frame named foo:
sink('some_output_file.txt')
uvals = unique(foo[,1])
for(uIDX in c(1:length(uvals))) {
    gname = uvals[uIDX]
    IDX = which(foo[,1] == gname)
    cat(sprintf("%s\t%s\n", gname, paste0(foo[IDX,2], collapse="\t")))
}
sink()

For that purpose, you may also use the reshape2 R package which allows you to switch wide to long format for data.frames. However, for that, you need an additional column with sample identifiers, e.g.:
G01000001_1 Sample1 0
G01000001_1 Sample2 0.00855717
G01000001_1 Sample3 0.0113467
G01000002_1 Sample1 0.0121086
G01000002_1 Sample2 0.0121108
G01000002_1 Sample3 0.0148402
G01000003_1 Sample1 0.0183057
G01000003_1 Sample2 0.0268394
G01000003_1 Sample3 0.0298547

Then, load and convert the data in R:
library(reshape2)
a<-read.table("test")
dcast(a, V1~V2, value.var = "V3")


text:
Error: FASTA file not found using building kallisto index
I am trying to build an index on kallisto on mac using R studio and have saved my FASTA file to my desktop, however, when I run kallisto index -i DI_mRNA.index DI_mRNA.fa.gz, there is an error saying Error: FASTA file not found DI_mRNA.fa.gz

In what way are you running Kallisto on R?
My guess would be you are not in the same directory as the fastq file. You can find out the current working directory with getwd() and set it with setwd().

As mentioned in a comment, your working directory lacks the fasta file.


text:
ERROR: (subscript) logical subscript too long
I am trying to remove a row of a dataframe in R however because the dataframe is so large when I use:
rawCounts <- rawCounts[rawCounts != geneID, ]

I get the error:
ERROR in rawCounts[rawCounts != geneID, ]: (subscript) logical subscript too long

Is there another way of removing a specific row in a large dataframe that is not so memory intensive?
Thanks in advance.
I have quite a bit of memory, memory.limit() returns:
> memory.limit()
[1] 16318

I am answering my own question here.
I realised I should have been using
rawCounts <- rawCounts[rownames(rawCounts) != geneID, ]

Problem solved.


text:
DESeqDataSetFromTximport all(lengths > 0) is not TRUE
I am suddenly running into an error when running the DESeqDataSetFromTximport
txi.rsem <- tximport(files, type = "rsem", txIn = FALSE, txOut = FALSE)
dds <- DESeqDataSetFromTximport(txi = txi.rsem,
                                colData = SampleFile,
                                design = ~ Compound)`

using counts and average transcript lengths from tximport
Error in DESeqDataSetFromTximport(txi = txi.rsem, colData = SampleFile,  :
  all(lengths > 0) is not TRUE

any suggestions what causes this error?

I found this same error as of today using this code:
pad.gene.data <-
  tximport(files = file.list.for.gene.analysis,
           type = "rsem",
           txIn = FALSE,
           txOut = FALSE)

dds.gene.level <-
  DESeqDataSetFromTximport(txi = pad.gene.data,
                           colData = pad.characteristic.data,
                           design = ~ condition)

I had done this same differential gene expression analysis using transcript counts aggregated to the gene level, on a mac about 5 weeks ago, but as of this past week switched to a windows and re-installed R 3.6.3 and whatever the current version of tximport and DESeq2 is.
I was able to replicate the results I had gotten on the mac 5 weeks ago removing all rows where all gene lengths are greater than 0:
pad.gene.data$abundance <-
  pad.gene.data$abundance[apply(pad.gene.data$length,
                                1,
                                function(row) all(row !=0 )),]

pad.gene.data$counts <-
  pad.gene.data$counts[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

pad.gene.data$length <-
  pad.gene.data$length[apply(pad.gene.data$length,
                             1,
                             function(row) all(row !=0 )),]

So first removing any gene where there is a length count of 0 from the abundance and counts matrix, then removing the genes from the length matrix so that they are all the same size. It was able to run that way.
I am definitely skeptical of why I would need to do that as it doesn't seem that it was advertised, if this is a new version, that we need to preprocess the data beforehand. I hope someone who might know more of the code can comment on this.

I had this issue as well. It looks like all zero-length genes also have zero expression across all samples. The code below makes sure only to exclude genes that have both zero length and zero expression.
txi = tximport(file_paths, type = "rsem", txIn = TRUE, txOut = TRUE)

zero_length_and_unexpressed = (apply(txi$abundance, 1, max) == 0) &
                              (apply(txi$length, 1, min) == 0)

txi$length = txi$length[!zero_length_and_unexpressed,]
txi$abundance = txi$abundance[!zero_length_and_unexpressed,]
txi$counts = txi$counts[!zero_length_and_unexpressed,]

dds = DESeqDataSetFromTximport(txi, sampleTable, ~ 1)
```


text:
Generating taxonomic hierarchy by species/genus name
I am attempting to create a reference library of DNA plant barcodes in the ITS2 barcode region for plants from a specific region (Panama). I downloaded all the sequences for plants that resulted in a search for “ITS2” “Panama” and “Viridiplantae” on NCBI GeneBank. For the purpose of assigning sequences of unknown plants (obtained from environmental DNA studies) based on similarity to the known sequences I want to clean up my data and present each sequence with its corresponding taxonomic ranks - Kingdom ; Phylum ; Class ; Order ; Family ; Genus; Species.
My problem is that the taxonomic information for each sequence upload in NCBI is all over the place –e.g. some sequences have additional taxonomic ranks listed in their information and there is no clear definition of what corresponds to which taxonomic rank. Since there are a few hundred species I want to refrain from filling in the information manually and I hoped to solve this problem by mining the taxonomic ranks I need from a different database according to the species or genus name in the sequence information. However, I have yet to find a dataset that is comprehensive enough to include all the plants and is simple to access.
Do you have any suggestions on how to achieve my goal?
Thanks.
According to a recent study, a mini DNA barcode region located on the ITS2 showed quite promising results for plant discrimination:
Moorhouse-Gann et al.(2018) New universal ITS2 primers for high-resolution herbivory analyses using DNA metabarcoding in both tropical and temperate zone. [link]


text:
Where can I find lncRNA expression data for different cell types?
Are there any publicly available databases providing expression data for long non-coding RNAs (lncRNAs) across cell types of multicellular organisms? Alternatively, are there lesser known UCSC tracks for this?

Example 1: I want to compare tissue-specific expression of the human lncRNA Xist across HeLa, GM12878, and HEK293 cell types [or
more].

Example 2: I want to compare tissue-specific expression of the plant lncRNA COLDAIR across cortex, phloem, and trichoblast cell
types [or more].

As long as the lncRNA is annotated you can use https://gtexportal.org/home/gene.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.

This is a fantastic resource for comparing lncRNA expression (or any annotated gene expression) across primary, non-diseased human tissues! Thanks for sharing, Mary.
Hopefully similar portals/databases will be developed in the near future for other multicellular organisms [and maybe even for cell lines?!].


text:
Removing stop codons from alignment file for dn/ds calculation
This is one of the question which is very similar to my problem or almost same to my problem
I ran with dummy data nucleotide alignment and the amino file it works perfectly fine,i can run the dn/ds from phangorn without any issue.
But when i ran with my data nucleotide alignment and the amino file I do get this warning
In kaks(carps) : sequence lengths are not a multiple of 3

and error.
The code
Any help or suggestion would be highly appreciated.

I would align your sequences with the software available here: https://github.com/veg/hyphy-analyses/tree/master/codon-msa

That is a good approach and very up to date answer.

thank you i will try this out as well

post-msa.bf  is it like shell script? im seeing this extension for the first time


text:
Combine structures in the same session - PyMOL
I'm a student currently looking at antibody responses against a viral target protein of interest.
I have my own, annotated PyMOL session of my protein and I also have .pdbs of crystallised antibody structures that bind to this protein and rougly know whereabouts it binds. I would like to know if it is possible to combine these two structures into a single session to observe or at least overlay the two structures together. Many thanks!

Is the question "Can I overlay/combine two structured in a single PyMOL session?"? (took me moment to figure what is actually the question, so I wanted to suggest how to make the question clearer).

Could you please describe what you did so far to solve the Problem?

Open the PyMOL session.
Use the load command to load your PDB file(s): load file.pdb.
Use the align or super commands to overlay structures, e.g. align prot1, prot2.


text:
What is a good RNA seq normalization method that allows for across sample comparisons and between transcripts
What is a good RNA seq normalization method that allows for across sample comparisons, and allows between transcripts comparisons as well? I read that TMM for example allows across sample comparisons but not between transcripts. On the other hand, I read that TPM is not good enough for between sample comparisons

I found a recent technique that allows for both, called geTMM.

I'm having a bit of trouble understanding the question, but we've used DESeq2 for doing both between-sample and between-transcript comparisons.
When comparing transcripts for short reads I use the variant-stabilised transformation, adjusted for transcript length. This is good for a visual overview of results (e.g. in a heatmap).
When comparing samples I use the variance-adjusted DESeq2 test, producing normalised $\textrm{log}_2FC$ values. This works well for quantitative results.
See this paper, where we use both of these approaches (e.g. Figure 2 shows both in a side-by-side view).

As mentioned by @gringer DESeq2 is a popular choice for normalization of counts between samples. However, because each gene is adjusted individually it is not well suited to compare within a sample.
For this a global normalization like TMM from edgeR would be more appropriate.
You can find more information in this tutorial from the hbc at Harvard.


text:
Understanding the "relativeness" in RNA seq experiments
Let's say we have two samples, sample 1 and sample 2, and we're running an RNA seq experiment. Consider genes A and B for both. For simplicity, let's say that when we run the RNA seq experiment, we get the following results:
gene A sample 1 count: 2
gene B sample 1 count: 4
gene A sample 2 count: 2000
gene B sample 2 count: 4000
Based on my understanding, RNA seq gives relative results. So I guess in this guess, how can you tell whether either one of the following are true:
1- You sequenced sample 2 a thousand fold more compared to sample 1
2- Sample 2 really has a thousand fold more expression in both genes A and B

We cannot be sure about this. For all we know, they could have been sequenced to the same depth and there is a biological difference that drives the 1000x (ridiculously large) fold change.
See above

The only thing you can say for sure is that gene B is expressed twice as much as gene A in samples 1 and 2. That's what "relative" means - relative expression of genes within a sample.

thanks this was helpful. just want to see what others have to say as well before I accept the answer

Oh of course. I’m no expert at this, but this is what I’ve encountered in our experiments.

You obviously have to normalize for differences in depth followed by a proper statistical test to get fold changes and significances. Sequencing depth is a technical confounder that has to be removed, and that is what normalization does. You have been posting a lot of very basic RNA-seq questions, I suggest you read the RNA-seq Hitchhikers guide: https://peerj.com/preprints/27283/ before posting new questions first.

The process of RNA-seq is a sampling process of all the RNA molecules. We cannot make reliable comparisons before we know how deep we sample. And since we know it is a sample from a population, the absolute value does not make sense.

So can we know how deep we sample? And knowing that, what would that allow us to do exactly? And how does this tie to sequencing "depth"?

This is usually the sequencing depth when I think about this problem.


text:
How many markers or loci are used in forensic DNA profiling and the effect of the DNA sequencing?
I am reading the Wikipedia page on DNA profiling, specifically STR analysis. It says that DNA profiling relies on matching 20 or so loci. It also gives an account of finding false identification. But that account was published in 2008.
What is current state of art on DNA profiling? How many loci are used? Do we still need a database for constructing the probabilistic distribution of matching?
I am particularly interested in the move toward DNA sequencing, how does that affect the profiling and matching probability calculation.
Suppose we sequence the whole available sample not restricted to particular loci.

When we have a single clean sample, is it true the probability of matching bounded by the probability of monozygotic twins resulting in the mismatching probability very close to $1/2$ of the monozygotic twin probability? The latter probability across all population is about $0.002$. It this upper bound reachable?

When there are admixture of several people/suspects' DNA in a sample and all the DNA (all 24 chromosomes) are intact.

When there are admixture of several people/suspects' DNA in a sample and some (the word covers the whole spectrum of degree/proportion) DNA's are damaged/missing.

In all the above cases, does the extraneous information such as gender, race/ethnicity affect the matching probability? I suppose the degree of effect increases in the order case 1 < case 2 < case 3.


text:
How to use the DI-tector software for defective genomes analysis?
I am trying to use a free .py script from this site
  http://www.di-tector.cyame.eu/

while unfortunately there's neither manual nor any documentation presented.
I have downloaded the .py script however I have absolutely no idea how it can be started. Can anyone give me an idea of how that can be done?
The corresponding article (https://pubmed.ncbi.nlm.nih.gov/30012569/) only displays parameters for executing it and there is no way to understand what steps must be performed firstly. Thank you.

No documentation is a major reason to not use this tool and to start looking for something else. The website seems to be in a poor state towards content, everything says "coming soon", I would be careful to build your analysis on this.

agree, while unfortunately, this is the only tool that I have found yet, that is dedicated to defective genomes analysis.

If you download the script and then type python DI-tector_06.py -h into the command line then you see explanations on the parameters. Maybe that helps getting started.

true, it says I have to do it this way: usage: DI-tector_06.py [-h] [-g HOST_REF] [-s MIN_SEGMENT] [-m MIN_MAPQ] [-n NB_READS] [-o OUTPUT_DIRECTORY] [-t TAG] [-d] [-l INDEL_LENGTH] [-f] [-p POLARITY] [-q] [-k] [-x NB_THREADS] Virus_Ref Input_Data Just cannot figure out what are Virus_Ref Input_Data might be just a .fastq and .gtf files?

Regardless it is a stylish (and fashionable) question

See the help:
python DI-tector_06.py -h
usage: DI-tector_06.py [-h] [-g HOST_REF] [-s MIN_SEGMENT] [-m MIN_MAPQ]
                       [-n NB_READS] [-o OUTPUT_DIRECTORY] [-t TAG] [-d]
                       [-l INDEL_LENGTH] [-f] [-p POLARITY] [-q] [-k]
                       [-x NB_THREADS]
                       Virus_Ref Input_Data

positional arguments:
  Virus_Ref             Virus genome reference sequence in FASTA format.
  Input_Data            File containing single reads in FASTQ format.

optional arguments:
  -h, --help            show this help message and exit
  -g HOST_REF, --Host_Ref HOST_REF
                        Host genome reference sequence in FASTA format.
  -s MIN_SEGMENT, --Min_Segment MIN_SEGMENT
                        Minimum segment length. Default is 15.
  -m MIN_MAPQ, --Min_MAPQ MIN_MAPQ
                        Skip alignments with MAPQ smaller than INT. Default is
                        25.
  -n NB_READS, --Nb_Reads NB_READS
                        Show only DVGs with counts reads > or egal to INT.
                        Default is 2.))
  -o OUTPUT_DIRECTORY, --Output_Directory OUTPUT_DIRECTORY
                        Enter a directory name that all compiled output files
                        will be saved in.
  -t TAG, --Tag TAG     Enter a tag name that will be appended before each
                        output file. Default is 'DI-tector'.
  -d, --DVG_sequences   Generate multi-fasta file with DVG sequences. Default
                        is (OFF).
  -l INDEL_LENGTH, --InDel_Length INDEL_LENGTH
                        Skip alignments with size of InDels smaller or egal to
                        INT. Default size is 1.
  -f, --Fasta           Select '-f' if data is in FASTA format fasta. Default
                        is FASTQ.
  -p POLARITY, --Polarity POLARITY
                        [0] Positive strand genome / [1] Negative strand
                        genome. Default is 0.
  -q, --No_Quantification
                        Inactive percentage quantification. Quantification
                        need bedtools Default is (ON).
  -k, --Keep_files      Keep intermediaite files (i.e. alignment etc...).
                        Default is (OFF).
  -x NB_THREADS, --Nb_threads NB_THREADS
                        Number of threads. Default is 1.


text:
about two group comparison of three group data in DESeq2 package
When we have three groups samples (A,B,C) with 25000 genes and the main interest is A vs B, should we limit samples only A and B for normalization to perform DEG analysis? Or better to include all samples (A,B, and C) for normalization and specify the comparison group within three groups?
res <- results(dds, contrast=c("tissu_type","A","B"))

No gold standard for this. If you probably are going to need C at some point then include it. If not then leave it out. If C is very different from A and B it might negatively influence normalization and dispersion estimation. Depends on your analysis goal. I would probably include it because if you ever need C you have to repeat the analysis and results from A/B would change.

I guess it makes more sense to include A and B. Imagine you are writing a paper, and it simply does not make sense if you include C, which you never talk about in the paper. But there is no absolutely correct answer for this. I also believe that the results are going to be very similar. If the results using two methods are so different, then there might be something wrong.

Thank you so much for your all comments!

As explained on the comments, having some other samples might distort your estimations, add an outlier or something alike.
The advantage of having all the three (or all) groups of samples before testing differences, because usually methods like limma, DESeq2 and others uses the information of variation of the genes across all samples to update their prior probabilities and the variation of the genes on the statistical model used to compare between groups.
If the method and your model are robust use all the three groups (generally the methods used are robust, but sometimes the models are not). By robust I mean that it doesn't miss any know (and available) source of variation. Otherwise the comparison might be comparing only the samples not the effect of the condition on the samples. If you only compare the samples it would be equivalent to use a t.test to the limma methods (not completely sure on the DESeq2 case). Or explained differently, you wouldn't take advantage of the statistical improvements implemented on these methods
Usually the reasons why these groups are sequenced on the same batch is because they are related (or to add some known ground truth), so you can release it with the data when communicating your results but you may omit that you used this third group on the methods section.

Plot out the PCA of your samples.  If, say, PC1 splits the C's away from the others, and PC1 is more than, say, 50% of the variance, or if the C samples are a lot more spread out than the A and B samples, keep it apart.  Otherwise, do them all together.
If the C's aren't too different from A and B, the math will benefit from including them.


text:
Non-exome variants called from whole-exome sequencing data
I'm working on rather standard whole-exome sequencing data and treat it the same as whole-genome sequencing data (aligning it to the full GRCh38 reference assembly and calling variants with no exome-related parameters).
I recently noticed that my variant calling results contain numerous variants mapped to a non-exomic region, incidentally, some of them were significantly associated with the phenotype I study.
My questions are as follows:

Are the non-exome variants in my exome data invalid? or could it be that the exome kit captured more than it was designed to capture?
Should I treat raw exome sequencing data any different than genome sequencing data, and if so, what parts of my pipeline should be modified to accommodate exome data?

These are called off-targets due to unspecific binding of DNA to the capturing beads during library preparation. This is normal and expected since capturing is not a perfect process. It is correct to align against the whole genome but you should simply limit variant calling to exons (or basically the on-target regions) and ignore everything else.

Thanks for the reply, makes perfect sense. Now, considering one of my significantly associated variant is from an off-target region, should I report it as a legitimate direction for further research? Consider that the variant has adequate coverage, passes strict filtering in the entire cohort and heavily implicated in my phenotype of interest.

I would also call variants a few bases out of the exons to cover splice site mutations.


text:
Grabbing all Immune related genes with databases in R
I am having trouble grabbing specific pathway info using databases in R. I have RNAseq results and I want to remove immune related genes from the current list that I have.
With a vector of gene names/symbols from RNAseq , I queried org.hs.eg.db to grab entrezIDs keys etc.
mapIds(org.Hs.eg.db,keys=mygenes,keytype ="ENTREZID",column="PATHID")

I thought about trying to use pathids from this database which I think are KEGG pathway ids which I can subsequently filter out the appropriate  ones but my list goes from 19,000 to 4,000 before I drop immune related genes due to a lot of NAs. I tried reactomedb, but the pathnames don't really let me filter by immune related or not.
Any advice is appreciated.

mapIds Bioconductor https://bioconductor.org/packages/devel/bioc/vignettes/AnnotationDbi/inst/doc/IntroToAnnotationPackages.pdf
KEGG https://www.genome.jp/kegg/
"KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies.
See Release notes (April 1, 2022) for new and updated features."

Hi James, those pathID come from the KEGG database and aren't updated since 2016 due to a license change. There isn't an easy solution to this problem (what is the threshold for a gene to be considered immune related?) and most people either use REACTOME, wikipathways or MSigDB (from the broad institute).

I was a bit confused on how to query REACTOME, as it was hard to visually see which terms were immune related when I looked at the pathname keys. Would a threshold use some sort of ranked pvalue (in which case probably something low)? Thanks for the suggestions on the others!

Well the problem is that the "olfactory pathway" on your epithelial cells might indicate an immune response and there isn't a test to prove that. People use some kind of network/system-biology approach or the databases available. Afaik there isn't much progress on this field. I'm attempting to tackle this problem by creating new methods about pathways but my thesis is in other topics so I can't invest much time on this...


text:
Finding out-group in unrooted tree
Im trying to find the out-group in the alignment my data files file a and  file b both the case i'm getting NA when i run this     pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]
I tried NJ instead of upgma i do still get NA not sure what is the issue ,but at the same time when i used other alignment file with 50 sequence  works fine.
a<-read.alignment("abhi_seq/final_abhi_tempa1.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]

the fact is that the present upgma tree as a midpoint rooting here, corresponding to node 46.
getRoot(upgma.a)

so for this kind of tree, I recommend you to use the midpoint() rooting function.
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
plot(midpoint(pars.a))

let me try it out and how do we decide which rooting to use? based on observation after plotting

yes visual inspection seems relevant. It totally depends of your input data, whether your sequences are strongly conserved (without evident outgroup) or more divergent. I think you may write a control loop after constructing the upgma tree, by checking if the ouptput of getRoot(upgma.a) returns a node or a tip. If it is a node, then it might proceed with the midpoint() version.

I haven't tried with ML but this mid point should work with ML as well?

yes, probably! it is probably a more universal solution, but not the necessarily the best one!

last time i was looking to compare tress which you had suggested turn them into dedrogram and then for comparison you suggested something this '
pars.a.rooted<-root(pars.a, outgroup="ABN10911", resolve.root=T)
pars.a.rooted.dd<-as.dendrogram(force.ultrametric(pars.a.rooted))' how to do the same  here shall i put node 46 as the outgroup?

It won't work because you have only 45 tip labels, if I remember well. I think you may run: pars.a.rooted.dd<-as.dendrogram(force.ultrametric(midpoint(pars.a)))

works perfectly fine for parsimony and...for pars for ML i m doing this rectify me if im wrong "ml.a.rooted.dd<-as.dendrogram(force.ultrametric(midpoint(fit$tree)))
"

seems ok. It is true that in the previous code, we used  root(midpoint(fit$tree), outgroup="AAA64460", resolve.root=T) to force using a given tip as root. Using force.ultrametric(midpoint(pars.a)) should work too indeed.

Let us continue this discussion in chat.


text:
differential analysis of chip-seq data
I have several sets of chip-seq data. I called the peaks using Macs2. I am pretty new to the field and I will appreciate any help. I wanted to annotate the peaks and see which peaks are shared between them and which ones are different. What is the best way to do this? which package should I use?
Many thanks

For sharing, you can test the proportion of overlapping peaks using bedtools intersect or the find.overlap function in GenomicRanges packages in R.
What makes more sense is perhaps to merge all the peaks together using bedtools merge, and then count the number of reads/fragments in each peak for each sample using featureCounts. The output of this is a matrix, which is very similar to gene expression results. Then you can do some differential tests on this data.

In a similar question, I recommended that directly comparing called peaks can be somewhat misleading. This topic is addressed in the documentation of the DiffBind R package on Bioconductor, in the "Comparison of occupancy and affinity based analyses" section.
But if you're just looking for genome annotation software, the annotatr package on Bioconductor is pretty useful

It depends on how many samples you have.
If you only compare one sample to the second, you can pool and call peaks and then directly calculate RPKM ratios for each peak.
It you have more samples, you may want to use IDR to score differential peaks.
You may also use packages such as Fseq or HOMER to treat one group as ChIP experiment while the other group as control to call peaks.
Alternatively, you could pool all the samples to call peaks, then calculate read counts in each peak for each sample, and then use DESeq2 or other differential tools to calculate significantly different peaks’ significance and effect size.


text:
output Array[Pair[X,Y]] to file inside a WDL task
I have a table containing different metadata information about a set of samples as an Array[Object] and I need to extract a couple of columns from this table and output those two columns into a file that will be read and used by the task. I wanted to know whether I can do this within a task rather than within a workflow, that is, without using the scatter construct which is not allowed inside a task.
If I zip two arrays I obtain an Array[Pair[String,String]] variable.
Now if I want to output to file an Array[Array[String]] I can use write_tsv(), if I want to output to file a Map[String,String] I can use write_map().
But WDL does not let me coerce Array[Pair[String,String]] to either Array[Array[String]] nor Map[String,String] (though I can coerce Map[String,String] to Array[Pair[String,String]]), so neither option is available.
Is there a solution to this that does not involve the scatter construct?

Hi @freeseek, could you please explain the biological relevance behind you informatics WDL question? Are you using BioWDL?

write_json() is supposed to provide the general way to serialize any data structure for processing/manipulation in the task command script (perhaps leveraging jq). Unfortunately, last I heard, Cromwell's implementation doesn't work on some useful cases right now. (Possibly out of date; I haven't tried lately)

Just learned this yesterday: in the development version you could try to use as_map() and then write_map(). However, you have to assume that all the left element in the Array[Pair[String,String]] array are unique or as_map() will fail.


text:
Scatter tasks around predefined batches inside a WDL workflow
Suppose I have two arrays Array[File] X and Array[String] Y both of length N, with Y being a sort of R-like factor that defines batch membership for X, for a total of M batches with M < N. A valid example would be:
Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

In which case N=5 and M=2. I want to scatter M tasks so that each task takes only the files corresponding to the same batch. In the previous example, one task would take as input:
Array[File] X =["file1", "file2", "file4"]
String Y = "A"

And the other task would take as input:
Array[File] X =["file3", "file5"]
String Y = "B"

Is there an easy way to do this within a WDL workflow in the general case? I understand that there are less elegant ways to do this as its own task but, as this seems like it could be a general purpose scatter/gather problem, an easy solution might be of general interest to the community.

Here's an alternative to using a custom task, albeit still rather unappetizing:
workflow groupby {
  input {
    File? nil
  }

  Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
  Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

  scatter (y in ["A", "B"]) {
    scatter (p in zip(X,Y)) {
      File? x_or_nil = if (p.right == y) then p.left else nil
    }
    Pair[String,Array[File]] group = (y, select_all(x_or_nil))
  }

  output {
    Array[Pair[String,Array[File]]] ans = group
  }
}

(It'll be slightly cleaner with the upcoming None literal instead of dummy nil)
FWIW, here's a lengthy discussion about the idea of allowing developers to define custom library functions which just run a small task that manipulates data structures (hinting that it needn't necessarily offload to a worker node, with the huge latency that typically entails).

When you use y in ["A", "B"] you are doing something that assumes you already know the unique contents of Y. This could be achieved using the keys() function but that is also only present in the development version and if the development version is to be used the solution with collect_by_key() is more minimalist.

It turns out the development version has a function collect_by_key that is particularly helpful for scattering around predefined batches:
Array[String] X  =  ["A",     "A",     "B",     "A",     "B"]
Array[File] Y = ["file1", "file2", "file3", "file4", "file5"]
scatter (pair in as_pairs(collect_by_key(zip(X,Y)))) {
  call test {input: X = pair.left, Y = pair.right}
}

This will call task t first with input:
String X = "A"
Array[File] Y =["file1", "file2", "file4"]

And then with input:
String X = "B"
Array[File] Y =["file3", "file5"]

This is a readable solution to this problem allowed by collect_by_key. This function is not very well advertised (I had completely missed it) but kudos to who introduced it in the language. The only drawback is that collect_by_key is not available in version 1.0.


text:
Sophisticated algorithms in biological systems
I am looking for examples of (known) algorithms/data structures in biological systems. My current understanding is that these don't exist - or haven't been found - beyond a relatively low level of complexity, which would seem like strong evidence to me that evolution isn't able to come up with something as complicated (and multi-step) as, e.g., the Schönhage-Strassen algorithm.
I heard of Bloom filters for recognizing smells in fruit flies (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6304992/).

It would be helpful to know more about what you mean by "algorithms" or "data structures". Organismal morphogenesis could be looked at as an incredibly complex set of linked algorithms- people call it a "developmental program".

There are many examples of biological systems doing computation, according to my understanding:

Regulatory networks encode decision trees of significant complexity: wikipedia.
Alan Turing did significant work on organismal morphogenesis showing how simple chemical patterns of inducers can lead to complex outcomes, "Turing patterns".
All of organismal development, cell division, etc. are tightly regulated procedures of gene regulation and morphogenesis.
The inspiration for neural networks is brains.


text:
Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design
I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in in vitro cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps.
> library(edgeR)
> load("KC_Raw.RData")
> y <- DGEList(counts = KCraw.data[,2:25])
> keep <- rowSums(cpm(y) > 10) >= 3
> targets <- data.frame(rownames=colnames(KCraw.data)[2:25] ,
+                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
+                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
+                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
> Group <- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep="."))
> targets <- cbind(targets,Group=Group)
> targets
   rownames Time FieldN CultureN      Group
1   KC1_H2L 2DIC     FH       CL FH.2DIC.CL
2   KC2_H2L 2DIC     FH       CL FH.2DIC.CL
3   KC3_H2L 2DIC     FH       CL FH.2DIC.CL
4   KC4_H2H 2DIC     FH       CH FH.2DIC.CH
5   KC5_H2H 2DIC     FH       CH FH.2DIC.CH
6   KC6_H2H 2DIC     FH       CH FH.2DIC.CH
7   KC7_L2L 2DIC     FL       CL FL.2DIC.CL
8   KC8_L2L 2DIC     FL       CL FL.2DIC.CL
9   KC9_L2L 2DIC     FL       CL FL.2DIC.CL
10 KC10_L2H 2DIC     FL       CH FL.2DIC.CH
11 KC11_L2H 2DIC     FL       CH FL.2DIC.CH
12 KC12_L2H 2DIC     FL       CH FL.2DIC.CH
13 KC13_H5L 5DIC     FH       CL FH.5DIC.CL
14 KC14_H5L 5DIC     FH       CL FH.5DIC.CL
15 KC15_H5L 5DIC     FH       CL FH.5DIC.CL
16 KC16_H5H 5DIC     FH       CH FH.5DIC.CH
17 KC17_H5H 5DIC     FH       CH FH.5DIC.CH
18 KC18_H5H 5DIC     FH       CH FH.5DIC.CH
19 KC19_L5L 5DIC     FL       CL FL.5DIC.CL
20 KC20_L5L 5DIC     FL       CL FL.5DIC.CL
21 KC21_L5L 5DIC     FL       CL FL.5DIC.CL
22 KC22_L5H 5DIC     FL       CH FL.5DIC.CH
23 KC23_L5H 5DIC     FL       CH FL.5DIC.CH
24 KC24_L5H 5DIC     FL       CH FL.5DIC.CH

I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example
> y <- DGEList(counts = KCraw.data[keep,2:25], group = Group)
> y <- calcNormFactors(y)
>
> TMM <- KCraw.data[keep,2:25]
> for (i in 1:24) {
+   TMM[,i] <- TMM[,i] / (y$samples$lib.size[i] * y$samples$norm.factors[i]) * 1e6
+ }
>
> y <- DGEList(counts = TMM,group = Group)
>
> design <- model.matrix(~0+Group)
> colnames(design) <- levels(Group)
> y <- calcNormFactors(y,method = "TMM")
> y <- estimateDisp(y,design)
> fitQL <- glmQLFit(y,design)
> fit <- glmFit(y,design)
> myKC.contrasts <- makeContrasts(
+   H2H.H2L = FH.2DIC.CH - FH.2DIC.CL,
+   L2H.L2L = FL.2DIC.CH - FL.2DIC.CL,
+   H2H.L2H = FH.2DIC.CH - FL.2DIC.CH,
+   H2L.L2L = FH.2DIC.CL - FL.2DIC.CL,
+   H5H.H5L = FH.5DIC.CH - FH.5DIC.CL,
+   L5H.L5L = FL.5DIC.CH - FL.5DIC.CL,
+   H5H.L5H = FH.5DIC.CH - FL.5DIC.CH,
+   H5L.L5L = FH.5DIC.CL - FL.5DIC.CL,
+   H2H.L2L = FH.2DIC.CH - FL.2DIC.CL,
+   H5H.L5L = FH.5DIC.CH - FL.5DIC.CL,
+   H5L.H2L = FH.5DIC.CL - FH.2DIC.CL,
+   H5H.H2H = FH.5DIC.CH - FH.2DIC.CH,
+   L5L.L2L = FL.5DIC.CL - FL.2DIC.CL,
+   L5H.L2H = FL.5DIC.CH - FL.2DIC.CH,
+   levels=design)
> design
   FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL
1           0          1          0          0          0          0          0          0
2           0          1          0          0          0          0          0          0
3           0          1          0          0          0          0          0          0
4           1          0          0          0          0          0          0          0
5           1          0          0          0          0          0          0          0
6           1          0          0          0          0          0          0          0
7           0          0          0          0          0          1          0          0
8           0          0          0          0          0          1          0          0
9           0          0          0          0          0          1          0          0
10          0          0          0          0          1          0          0          0
11          0          0          0          0          1          0          0          0
12          0          0          0          0          1          0          0          0
13          0          0          0          1          0          0          0          0
14          0          0          0          1          0          0          0          0
15          0          0          0          1          0          0          0          0
16          0          0          1          0          0          0          0          0
17          0          0          1          0          0          0          0          0
18          0          0          1          0          0          0          0          0
19          0          0          0          0          0          0          0          1
20          0          0          0          0          0          0          0          1
21          0          0          0          0          0          0          0          1
22          0          0          0          0          0          0          1          0
23          0          0          0          0          0          0          1          0
24          0          0          0          0          0          0          1          0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"

> myKC.contrasts
            Contrasts
Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
  FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
  FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
  FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
  FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
  FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
  FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
  FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
  FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
            Contrasts
Levels       L5H.L2H
  FH.2DIC.CH       0
  FH.2DIC.CL       0
  FH.5DIC.CH       0
  FH.5DIC.CL       0
  FL.2DIC.CH      -1
  FL.2DIC.CL       0
  FL.5DIC.CH       1
  FL.5DIC.CL       0

After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.
> myKC.contrasts <- cbind(myKC.contrasts,
+                         Development = c(1,1,-1,-1,1,1,-1,-1),
+                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
+                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
+ )
> myKC.contrasts
           H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
           L5H.L2H Development FieldN CultureN
FH.2DIC.CH       0           1      1        1
FH.2DIC.CL       0           1      1       -1
FH.5DIC.CH       0          -1      1        1
FH.5DIC.CL       0          -1      1       -1
FL.2DIC.CH      -1           1     -1        1
FL.2DIC.CL       0           1     -1       -1
FL.5DIC.CH       1          -1     -1        1
FL.5DIC.CL       0          -1     -1       -1

Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast.
> lrt <- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
> topTags(lrt,n=nrow(y$counts))["GRMZM2G445575",]
Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL
                 logFC   logCPM        F       PValue          FDR
GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
# FC is a data frame of the logFC of each constrast in columns for each gene in rows
> sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
[1] -6.636197

My first question is if this analysis is a valid way of summarizing the simple effects of each treatment.
I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors.

I don't quite get why you need to manually change the counts with a for loop. It is used as a size factor in the linear model. Now, for the pairwise comparison, what you have in the first part is ok

For the general culture N etc, run a linear model and run the estimateGLM... etc stuff, it should give you the effects you want. You cannot get everything from one model

I think the problem is in the design. There is no room for error or variation from a common ground.
I'm not sure of the solution, but I think you'll need to delete a column to give the design a degree of freedom, and probably add the intercept. You'll need to adjust the contrasts accordingly.
However, I would recommend to ask this at support.bioconductor.org. There are more experts on linear modelling and contrasts there than me. (If you ask there make it easier to copy and paste your code)

Thanks for the replies. I'm not really sure how to do the modifications you recommended to my design or contrasts. I am writing a question on bioconductor, but how can I make it easier to copy and paste my code? This was my first post on a bioinformatics support forum.

Ideally you could use a package as reprex. But if you use code blocks and remove the > at the beginning of the code, comment the output (so that if I paste it on the console it doesn't error) you'll do good. Another thing to consider is that I cannot load("KC_Raw.RData") because I don't have it, so try to make the example self contained. But the code provided is enough to understand your problem so don't worry if you use the same code used here there.


text:
How to make work programs from the $PATH?
I am trying to analyze my RNAseq reads for defective genomes and I use this program (http://www.di-tector.cyame.eu/) that is suggested by Beauclair et al (https://pubmed.ncbi.nlm.nih.gov/30012569/)  for that.
After downloading a .py script and running it with the following command:
sudo python3 DI-tector_06.py /mnt/e/nastya/SLX066-02/sequence.fasta /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz

It says at the beginning that
=================================
Program: DI-tector
Version: 0.6 - Last modified 25/05/2018
=================================
Requirement: (must be in your $PATH)
-bwa
-samtools

Optional: (must be in your $PATH)
-bedtools

I have tried to install bwa using:
sudo apt-get update -y
sudo apt-get install -y bwa

and samtools using
cd /usr/local/bin
sudo wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
sudo tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
sudo make

I did not get exactly what does that mean: (must be in your $PATH) but I have searched here and there and discovered that it might mean /usr/local/bin
The program starts and, after some calculations it gives me the following error:
Input file:               /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz
Host reference:           None
Virus reference:          /mnt/e/nastya/SLX066-02/sequence.fasta
Remove segment:         < 15 nt
Remove reads with MAPQ: < 25
Allow InDel. length:    > 1 nt

=================================
Reading gziped input file...
=================================
Number of reads in input file: 18,035,544
=================================

=================================
Step 1/5 : Alignment against Viral reference...Filtering...
=================================
[E::bwa_idx_load_from_disk] fail to locate the index files
Traceback (most recent call last):
  File "DI-tector_06.py", line 183, in <module>
    Popen(["samtools", "view", "-Sb", "-F", "4"], stdin=out_file_onVirus, stdout=out_file_Virus).communicate()[0]
  File "/usr/lib/python3.6/subprocess.py", line 709, in __init__
    restore_signals, start_new_session)
  File "/usr/lib/python3.6/subprocess.py", line 1344, in _execute_child
    raise child_exception_type(errno_num, err_msg, err_filename)
PermissionError: [Errno 13] Permission denied: 'samtools'

I have no idea how samtools must be accessed, I tried sudo command with python3 but it did not work.

generally, if you type sudo make install, after you run sudo make, it will install the binary in your $PATH and you won't need to worry about adding it there afterwards

A piece of advice: Dont "over-sudo": You shouldn't need sudo to run a data analysis program. When installing samtools, avoid downloading stuff in /usr/local/bin: Do it in a local directory (I would recommend having a $HOME/src for that), and use sudo only to install the program to /usr/local/bin.

"Permission denied: 'samtools'" may mean that the samtools "executable" is actually not executable. This can be because it is on a disk that doesn't allow execution of files, but more likely it just needs a fix of the file permission (chmod +x <path_to_samtools>).

echo $PATH will show the directories in your search path. Either copy the samtools executable there or add its location to the $PATH variable. You would be well served with a unix command line course or tutorial.

you can add a directory to $PATH by executing following command in the terminal:
export PATH="$HOME/Documents/Software/bbmap:$PATH"
In this example, all executable commands in the folder bbmap can be executed from any other directory in the terminal.
So you can install samtools in an appropriate directory and then execute export PATH="$HOME/your/directory:$PATH". Then check if it succeeded by running echo $PATH. This latter command will return a list of directories which are assigned to $PATH and your new directory (where samtools is located) should be included.

You can run a piece of software by specifying exactly where it is, like /home/usr/software/samtools_1.1.0/samtools
Or if that location is in your $PATH, you can run it with just samtools.  The software you downloaded needs to run those programs, but there's no way for you to give it the specific path for them; it's going to try and run them with just samtools.  So you need to make sure that will work, by putting their locations in your PATH.
It is not clear to me that running every command with sudo is wise.  If you don't understand what $PATH is, I'm pretty sure you should not be running everything with sudo.

The clean way to do this is to make your own bin within your user directory
cd ~
mkdir bin
vi .profile # or vi .profile_bash or for OS X using `csh` or `tcsh` vi .cshrc

For bash (Linux, or OS X) within vi type "i" then type,
PATH=$PATH:~/bin

or
PATH=~/bin:$PATH

... followed by esc key and "wp" enter
Then enter for csh or tcsh e.g. for OS X
set path = ($path /Users/username/bin /Users/username/anaconda3/bin)

Same instructions for using vi as for bash. Dump bedtools into any directory (not the desktop) and connect the programs into the bin using ln
cd ~/progs/bedtools
ln -si $PWD/prog ~/bin

If you are using csh or tcsh you need to type rehash but bash does this automatically.
The programs will then be cleaning connect to your $PATH via ln without massive memory usage and visible anywhere within your home space. Its good practice for any Linux / OS X (time gone by UNIX) program or script and whilst you don't need it for this operation as you add more stuff its better.


text:
.SDF file to .PDB file
I am currently working on identifying effective antiviral compounds against the corona virus. I have downloaded the Covid-19 AntiViral Candidate Compounds Dataset from the www.cas.org webpage, the link to which I have provided below. Me and a few other doctors are using Docking software to predict the docking capacity of these anti-viral compounds. Docking software requires the anti-viral compounds to be in the .pdb file format only then is it recognized by Automated docking software. Is there any way I can convert the .sdf file containing the anti-viral compounds to .pdb file or is there a way of obtaining the .pdb files for list of anti-viral compounds.
Thank you for your time and patience,
Anas Hussain.
Anti-Viral Candidate Compounds Dataset
https://www.cas.org/covid-19-antiviral-compounds-dataset

Use OpenBabel.
Something like babel -isdf file_in.sdf -opdb file_out.pdb.

Bizarrely, Obabel seems to have written a PDB that was incorrectly spaced for the OP as seen in Error pasing the following line in pdb

If you can use python, rdkit is a great package for cheminformatics —best installed via conda though.
It allows you do all sorts of things, such as match substructures, embed conformers, edit compounds, compute advanced parameters and many forms of shape & colour finger-/foot-printings.
Given a multi-compound sdf you can read each individual molecule and write them out as PDBs.
from rdkit import Chem
with i, mol in enumerate(Chem.SDMolSupplier('in.sdf')):
    Chem.MolToPDBFile(f'out_{i}.pdb')

We can be fancier with the output name
from rdkit import Chem
import re
with mol in Chem.SDMolSupplier('in.sdf'):
    # The headerline of a mol file is retrieved oddly.
    name = mol.GetProp('_Name')
    # make filename safe
    slugify = lambda name: re.sub(r'[\W_.-]+', '-', name)
    Chem.MolToPDBFile(f'{slugify(name)}.pdb')

To specify what name (e.g. CA) each atom gets and what is the residue name, chain and residue position —information absent in the sdf, it gets a bit odd.
  names = [' CA ', ' CB ', ' CG ', ' CD '] # say
  for i, atom in enumerate(mol.GetAtoms()):
        info = Chem.AtomPDBResidueInfo(names[i], # or whatever.
                                       residueName='LIG')

For more see this blog post.

for the rdkit sdf to pdb conversion, should i enter the multi entry sdf filename in the line below;
"with mol in Chem.SDMolSupplier('in.sdf'):"

Correct in.sdf is an example filename —the quotes must be there or it will think it a variable. If you are unfamiliar with python this route is probably not the best way.

I am unfamiliar with python however I am currently trying to learn it from scratch.
Thank you


text:
Plotting metadata from GenomicRanges object
I have a GRanges object with metadata that I want to plot. Many of the online tutorials which manipulate these GRanges for plotting show how to plot the main columns (chromosome,start/end,strand) but nothing to plot the metadata as a function of the position in each row.
How precisely should I do this? Ideally, I would like to plot metadata from multiple GRanges objects on the same plot to compare them.
EDIT (following comment):
GRanges object with 6 ranges and 140 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1 100000-199999      * |
  [2]     chr1 400000-499999      * |
  [3]     chr1 600000-699999      * |
  [4]     chr1 700000-799999      * |
  [5]     chr1 800000-899999      * |
  [6]     chr1 900000-999999      * |
                  eigen   domain         arm
              <numeric> <factor> <character>
  [1] -1.62475825405065     open          1p
  [2] -1.62475825405065     open          1p
  [3] -1.62475825405065     open          1p
  [4] -1.62475825405065     open          1p
  [5] -1.63354915136132     open          1p
  [6]  -1.6322733358426     open          1p
      exon.fraction       gc     short      long
          <numeric> <matrix> <numeric> <numeric>
  [1]       0.02968  0.45815         0         0
  [2]       0.04305  0.40192         0         0
  [3]       0.00959  0.42091         0         0
  [4]       0.04593  0.43337         0         1
  [5]       0.02425  0.44822         2         1
  [6]       0.16977  0.62123        27        56
         nfrags             ratio   short.corrected
      <numeric>         <numeric>         <numeric>
  [1]         0               NaN              <NA>
  [2]         0               NaN              <NA>
  [3]         0               NaN              <NA>
  [4]         1                 0 -3.50822623240367
  [5]         3                 2  2.81796238367926
  [6]        83 0.482142857142857  39.8824095079959
...

My question regards plotting a column like eigen as a function of the position (given by seqnames and ranges)

Can you show us head(., 10) for a few of the datasets so we can understand what you're trying to plot?

You can access the data.frame just using $. All of them can be accessed using values() . Something like this:
library(GenomicRanges)
gr = GRanges(seqnames=1,IRanges(start=1:10,width=1))
values(gr) = data.frame(eigen=rnorm(10),ratio=runif(10))

plot(start(gr),gr$eigen)

or
plot(start(gr),values(gr)$eigen)


text:
Do I need to remove water molecules before protein-ligand docking?
I am going to screen small molecules against a protein target. Do I need to remove all the solvent molecules before performing VLS?

Yes, if you are using implicit solvent you ought to remove all of them.
Exception. Unless the crystallographic waters are trapped essential structural components of your active site and your ligand is not going to replace it, in which case keep it because its remove changes the energy of the active site —say you have a catalytic triad, with a water that is pinned between the acid (e.g. aspartate) and the base histidine, its removal would change the tautomerisation of latter and result in the nucleophile being protonated (and not nucleophilic). Having said that, if you don't have time to look at your structure (or cannot because of automation) and are okay with some artefacts from incorrect protonation strip them.
If you are not using implicit solvent, then it depends on the tool used —most require them stripped and re-added. If you do keep them, make sure that the crystallographic water (HOH) has the same residue code as the added water (TP3, WAT etc.).


text:
getGeneLengthAndGCContent in EDASeq complains about "Ensemble gene ID"
I am trying to use the getGeneLengthAndGCContent function from the EDASeq library to retrieve gene lengths for c. elegans. I retrieved Ensemble gene ids from Biomart for c.elegans and I'm using those as input to getGeneLengthAndGCContent (I picked "ensemble_gene_id" from the filters of the worm annotation Mart). Here's an example of how I'm trying to use the function:
getGeneLengthAndGCContent("WBGene00001042", org, mode=c("biomart", "org.db"))

When I try to run this, I get
Error in getGeneLengthAndGCContent("WBGene00001042", org, mode = c("biomart",  :
  Only ENTREZ or ENSEMBL gene IDs are supported.

@ATpoint oh that's weird why does BioMart list it as such?

@ATpoint also when I look at Ensembl I don't find any other identifier for a gene of c elegans (here's an example http://uswest.ensembl.org/Caenorhabditis_elegans/Gene/Summary?g=WBGene00001042;r=III:5415556-5419565;t=W03A5.7.1). Could you guide me as to where to find those other ids?

Never mind, I was wrong.

@ATpoint so WBGene00001042 is indeed the Ensembl id?

I found my answer here. It uses getGeneLengthAndGCContent but changes a small part of the function.


text:
How to calculate module-trait relationship when trait data is in binary format?
I have a dataset of 50 breast cancer samples. These samples are classified into four subtypes Lum A, Lum B, Her2 and Basal. I have been working with lncRNAs and protein-coding genes. To identify the functions of lncRNAs, I have used WGCNA through which I get the modules with protein-coding genes and lncRNAs. The expression data is stored in a matrix datExpr.
I have the trait information in datTraits like below:

I totally got 20 modules named with different colors stored in dynamicColors. Like below I first calculated eigenvalues and then used it and trait information to check the association between Modules and trait data.
MEs0 = moduleEigengenes(datExpr, dynamicColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);

What is the optimal method to associate modules and trait using binary data?

Is Pearson correlation a right way to check the association?

Is there any other way to check the association between modules and trait (which represented like above in binary)? Details of methology would be welcome.

Patients are not genetically related and I have no information on environmental risk factors. I was thinking about logistic regression, but before doing that, I would like to take suggestions.

Answer from @m, converted from comment:
Personally I'd build a tree. Whats your question patients or the relationship between RNAseq 'signatures'? Its not the definitive analysis, but its good mining.
Pearson is doable but how many Pearson's? You'll end up with a 2x2 matrix. How will you resolve the matrix clustering? A maximum likelihood phylogeny solution would work here making use of the binary format. There are other approaches.
Spearman correlation would be better than Pearson's you end up with a 2x2 matrix.
Answer from @peter-langfelder, converted from comment:
You can use Pearson correlation, it is equivalent to carrying out an equal-variance t-test for the module eigengene between the two groups specified by the binary variable.


text:
classifying samples by TCGA signature
I have some RNA-seq samples from multiple glioblastoma tumours that I'm now trying to classify according to a specific gene signature (from Verhaak et al., 2010) using R. The gene signature is reported as a gene list with specific centroids for each of the 4 clusters (https://api.gdc.cancer.gov/data/941f81a1-05d7-4f84-80ec-534b8dc1ebac). I'm wondering how I can use this signature to classify my samples in R? Would it involve some sort of k-nearest neighbours method?
Additionally, the signature was identified using microarray data, but I am classifying RNA-seq data. Is there any sort of adjustment I should make to the signature to account for this?
Thanks in advance!

Hi inteeresting post and right site to post it :-)  Coudl you summarise the genetic signature here rather than use a weblink please?

Hi Michael! Thanks for the response. The authors of the paper analyzed 200 glioblastoma samples, and used consensus average linkage hierarchical clustering to identify 4 clusters in their dataset. Next, they used ClaNC, a nearest centroid-based classifier that balances the number of genes per class, and identified signature genes for all four subtypes. The gene signature is a collection of 840 genes (210 genes per subtype). I'm not familiar with gene signatures, but for each gene they have 4 centroid values corresponding to each of the subtypes.

This is not a very simple problem. I saved the table and "test.txt" and visualize it below:
library(pheatmap)
x = read.delim("test.txt")

head(x)
  Gene.Symbol            Composite_chr_coords  Proneural      Neural
1      CDKN1B   [36.1:12:12761576-12766570:+] 0.03353670  0.03353670
2       EPB41    [36.1:1:29086215-29319545:+] 0.03882308  0.03882308
3        CLGN  [36.1:4:141529057-141568232:-] 2.12796056  0.41999847
4      PDE10A  [36.1:6:165664525-165995575:-] 1.66012916  0.25349464
5     RALGPS2  [36.1:1:176960923-177152465:+] 0.78971500 -0.01945154
6        TAF5 [36.1:10:105117714-105138812:+] 1.08354870  0.09461930
    Classical Mesenchymal
1  0.03353670 -0.84379400
2 -0.84821990  0.03882308
3  0.41999847  0.41999847
4  0.25349464  0.25349464
5 -0.01945154 -0.01945154
6  0.09461930  0.09461930

You can see there are no so called markers genes and if you look at the scale of the data, it was normalized in microarray and it will be very hard for you to get your RNA-seq to fall into that range, and apply the medoids to identify your clusters or type.
Since you did not provide any data, and it is not trivial to simulate data for this, I can suggest the approach in words here:

Subset your count table to these exact set of 840 genes.

Normalize the data and do a PCA or hierarchical clustering to see if you indeed observe 4 or 3 or 2 subtypes

If you indeed see 4 or 3 clusters this would be straight forward, take the centres of means of these clusters, and correlate with each of the sub types here. Chances are the subtype that is most correlated with one cluster, is the subtype you are interested in.