text:
Annotation label not matching to the color assigned in complex heatmap
The issue im facing is I'm assignig the colour to the annotation label but when it is rendered the label in my metadata does't match with the final output not sure what exactly is the problem.
The code
data file  ka_ks and metadata
As an example this sample Q6Q1S2 in my metadata is Alpha but in the figure it is labelled under Gamma
Any suggestion or help would be really appreciated
    data <- read.csv("ka_ks.txt",header = TRUE,row.names = 1,sep = ",")
head(data)

#metadata <- read.csv("NUC_META_AMINO_heatmap_metadata_sans_nucleotide.txt",sep = ",")
metadata <- read.csv('vir_meta.txt',sep = ",")
#head(data)
head(metadata)
#names(metadata)[3] = "Host"

#Set annotation
ann <- metadata[,c(2,3)]
colnames(ann) <- c("Class","Host")
#rownames(ann) = metadata$Sample
ann
#ann = metadata
colours <- list("Host"=c("Human SARS coronavirus 2"="green4",
                         "Human SARS coronavirus"="darkgreen",
                         "Bat coronavirus"="wheat4",
                         "Murine coronavirus"="turquoise","MERS"="orange",
                         "Bovine coronavirus"="coral",
                         "Rat coronavirus"="saddle brown",
                         "Human coronavirus"="dark green",
                         "Porcine hemagglutinating encephalomyelitis virus"="maroon",
                         "Avian infectious bronchitis virus"="dark slate gray",
                         "Porcine respiratory coronavirus"="slategrey",
                         "Canine coronavirus"="purple",
                         "Porcine transmissible gastroenteritis coronavirus"="pink",
                         "Feline coronavirus"="gold",
                         "Porcine epidemic diarrhea virus"="brown"),
                "Class"= c("Alpha" = "blue" ,"Beta" ="black","Gamma" ="red"))

heat <- scale((data))
head(heat)

colAnn <- HeatmapAnnotation(df=ann, which="col",col=colours, annotation_width=unit(c(1, 4), "cm"), gap=unit(1, "mm"))
myCol1 <- colorRampPalette(c("navyblue", "white", "red"))(100)
myBreaks1 <- seq(-4, 4, length.out=100)

#Add column boxplot
boxplotCol <- HeatmapAnnotation(
  boxplot = anno_boxplot(data.matrix(heat),
                         border = TRUE,
                         gp = gpar(fill = '#CCCCCC'),
                         pch = '.',
                         size = unit(2, 'mm'),
                         axis = TRUE,
                         axis_side = 'left',
                         axis_gp = gpar(fontsize = 12)),
  annotation_width = unit(c(1, 5.0), 'cm'),
  which = 'col')

hmap1 <- Heatmap(heat,
                 name = 'ks_score',
                 #name="Z-score",
                 col=colorRamp2(myBreaks1, myCol1),
                 heatmap_legend_param=list(color_bar="continuous",
                                           legend_direction="vertical", legend_width=unit(5,"cm"),
                                           title_position="lefttop", title_gp=gpar(fontsize=10, fontface="bold")),

                 #Split heatmap rows by gene family
                 #split=df$Family,
                 #split = NULL,

                 #Row annotation configurations
                 cluster_rows=TRUE,
                 show_row_dend=TRUE,
                 row_dend_side = "left",
                 #row_title="Transcript", #overridden by 'split' it seems
                 row_title_side="left",
                 row_title_gp=gpar(fontsize=30,fontface="bold"),
                 row_names_gp = gpar(fontsize = 10),
                 show_row_names=TRUE,
                 row_names_side="right",
                 row_title_rot=0,

                 column_dend_reorder = TRUE,
                 #row_dend_reorder = TRUE,

                 #Column annotation configuratiions
                 cluster_columns=TRUE,
                 show_column_dend=FALSE,
                 column_title="Samples",
                 column_title_side="top",
                 column_title_gp=gpar(fontsize=20, fontface="bold"),
                 column_title_rot=0,
                 column_names_gp = gpar(fontsize = 10, fontface="bold"),
                 show_column_names=TRUE,

                 #Dendrogram configurations: columns
                 #clustering_distance_columns="euclidean",
                 clustering_method_columns="ward.D2",
                 column_dend_height=unit(10,"mm"),

                 #Dendrogram configurations: rows
                 #clustering_distance_rows="euclidean",
                 clustering_method_rows="ward.D2",
                 row_dend_width=unit(15,"mm"),

                 #top_annotation_height=unit(0.5,"cm"),
                 #top_annotation=ColAnn,
                 top_annotation=colAnn)

#bottom_annotation_height=unit(2, "cm"),
#bottom_annotation=boxplotCol)

draw(hmap1 ,heatmap_legend_side ='right',annotation_legend_side = 'right')

fixed it ..the orders between my metadata and the data matrix was not matching


text:
How to extract unmatched reads using bwa and samtools?
I have a single read (NOT paired) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
Here in Materials and methods the procedure is described as the following:
The first step of the workflow consists of an alignment of the reads against the host genome (Fig. 2i).
This step aims at discarding reads that map to the host genome and may partially map to the viral genome after segmentation, and at reducing the working file size.
For example, MV and rMV-ΔV data sets were generated from total RNA samples of infected cells and contained mostly reads mapping the human genome (≈99%).
This step uses a combination of bwa mem and samtools view with the parameters –bS –f4.
An additional step consists of an alignment of the reads against the viral genome of interest, in order to exclude reads perfectly mapping the viral genome.
Therefore, only unmapped reads are further analyzed.
Of note, clipped reads (i.e., CIGAR motif contain S or H) are also conserved. Some of these reads may map to viral genome recombination junctions that are present in DI genomes.
I have already been suggested to use bowtie2 instead of bwa but, firstly, the output is not clear for me and secondly, I would like to test the official protocol.
Since the article suggests to use bwa and samtools on this very first step, that is what I have done so far:

(Optional) not sure if this is important but as suggested by someone I have transformed .fna --> .fa
cp GCF_000001405.39_GRCh38.p13_genomic.fna GCF_000001405.39_GRCh38.p13_genomic.fa

Indexed human genome with bwa
bwa index GCF_000001405.39_GRCh38.p13_genomic.fa

aligned my only read with indexed genome from step 1
bwa mem GCF_000001405.39_GRCh38.p13_genomic.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz > whole.sam

Transform .sam --> .bam
samtools view -S -b whole.sam -o whole.bam

Separated unmapped reads (as it is recommended in Materials and Methods using -f4)
samtools view -f4 whole.bam > sample.unmapped.sam

Converted unmapped reads into .fastq format (since this is the format used by the software later)
samtools fastq sample.unmapped.sam > unmatched.fastq

I wonder if these steps are correct?

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to:
bwa mem (...options) | samtools view -o out.bam
samtools fastq -f 4 out.bam > unmatched.fastq

Thank you! What about step 0 - converting from .fna to .fa?

Unnecessary since this is only a suffix, the content of the file is not changed. You could have simply renamed the file with mv. I am 99.99% sure bwa does not care about the suffix of the fasta.

Thanks! I also wonder what are these 5 additional files after indexing, e.g. *.fa.amb, *.fa.ann, *.fa.bwt, *.fa.pac and *.fa.sa ?

Internal files that bwa needs for alignment (suffix arrays, burrows-wheeler transformed genomes...), nothing that you as a user have to worry about. Most are binary anyway.

Hi @DmitriiTrubetskoy please remember to upvote ""ATpoint" for their answer and mark the answer as "accepted" if this correctly answers your question.

Thanks!!! .. :-)

source: Stack Exchange

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

sort: votes


text:
Generating 3D coordinates error
I have been trying to convert a multi entry sdf file to pdb using openbabel. Everytime I check the Generate 3D coordinates box I end up with .ENT files instead of .PDB files. What am I doing wrong.
Thank you.

Voting to close, lacks details and reproducibility.

An ent file is a pdb file. Just change the extension. PDBe for example provides their PDB files as such.

ok. I tried opening the .ent file with Autodock but I kept getting the error, .ent file is an unrecognized format. How can I simultaneously change the extension for 50,000 ent files to pdb without having to do it manually for each one individually.
Thank you,
Dr. Hussain.

In your file manager, select all and rename

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


text:
RNA-Seq data transformation prior to sample correlation analysis
If I’m starting with Deseq2 normalized counts what are some preprocessing steps that I should apply to these data before estimating sample correlation using the cor function in R? For example, would it make sense to quantile normalize the data?

Hi @Alexis your question doesn't quite make sense, you are starting with normalised counts and you want to further normalise them? If you data is normalised you can proceed with Spearman's correlation analysis.

You don't absolutely need to perform any further steps before calculating the correlation. However you would likely benefit from the following sort of workflow:

Transform your data with vst() or robust log (see the DESeq2 vignette) so genes are on a more or less similar scale.
Select the top 300-500 most variable genes
Use this subset with to compute sample-wise correlation.

This also what follows what you would find in the DESeq2 user guide, which I strongly encourage you to read.


text:
How do I set a neural network to loop multiple times and average the resulting values?
I have a script in R/RStudio which creates random datasets of binomial variables, feeds them through a neural network, and calculates their likelihood ratio statistic and deviance. I'd like the script to loop with the seed incremented by 1, and collect the value from each run into a data frame, from which the values get averaged. For example, if the script returns 5, 10, and 15 on the first, second, and third iterations of the loop, I'd like the script to create a data frame with those values and then calculate their mean.
Here is my script, annotated:
library(neuralnet)

x=1021

#change nrow to number of times code is looped
result <- data.frame(matrix(nrow = 1000, ncol = 2))
colnames(result) <- c("LogL", "D")

repeat{

set.seed(x)

#2 inputs, max number of terms 3
x1=rbinom(2000,1,0.5)
x2=rbinom(2000,1,0.5)
y=rbinom(2000,1,0.5)

data=data.frame(x1,x2,y)

head(data)

#hidden units = 2
nnet=neuralnet(y~x1+x2,data,hidden=2,err.fct = "ce",linear.output = F,likelihood = T)

predictions=compute(nnet,data[,1:2])

predictions$net.result

#likelihood ratio stat
argument=y*log(predictions$net.result)+(1-y)*log(1-predictions$net.result)

head(argument)

#sum of likelihood ratio stat
LogL=sum(argument)

LogL

D=-2*LogL

D

result[x-1022, 1] <- LogL
result[x-1022, 2] <- D

result

#results <- data.frame(LogL, D)

x=x+1
if (x==2021) {
  break
}

print(x)

}

result

My repeat loop works, but a thousand iterations, for example, takes multiple minutes. Is there a more efficient way to code such a loop?

Your repeat termination conditions will never be satisfied. Every round, x will be set to 1021, incremented by 1 and compared to 1023, which will never be true. How will the repeat ever stop?

There are other approaches to parameterisation in ANN. @RamRS question on terminating recursion is definitely a good question.

Thanks, edited that and it works now, but it's super slow. I'll update the question!

The concern with your approach is called 'leakage' because you are parameterising the same data set that you are training. This can easily lead to overtraining is in ANN (artificial neural networks) and is not cool. Overtraining is where the ANN very tightly wraps around the training data set that it performs poorly on any other data set. In other words it is overtightened.
The way they get around it is to split the data set usually 60:20:20 into a training split, parameterisation (validation) split and test (testing) split. If the parameterisation split gets mixed with the training set - leakage happens and overtraining results. Normally the data is trained, and the trained model parameters are assessed/optimised in the parameter split. The trained data and optimised parameters are then tested in the test split.
As long as you are aware of the approach you are using, that I guess that is ok. If you separate your training and parameterisation (validation) data that would be cool.


text:
How to design DESeq2 LRT model with individuals nested in 2 levels?
We have a complicated experimental design that we would like to perform LRT analysis for. Our main goal is to discover significant genes for the "Injection:Social" interaction term across the entire dataset by removing it from the LRT reduced model, and as a bonus we are also interested in discovering significant genes for that interaction term for each respective brain region.
Sample  Injection   Social  Region  Individual  ind.n
HY06    L   ISO HY  S06 S1
NST6    L   ISO NS  S06 S1
TN06    L   ISO TN  S06 S1
HY08    L   ISO HY  S08 S2
NST8    L   ISO NS  S08 S2
TN08    L   ISO TN  S08 S2
HY30    L   KF  HY  S30 S1
NST30   L   KF  NS  S30 S1
TN30    L   KF  TN  S30 S1
HY32    L   KF  HY  S32 S2
NST32   L   KF  NS  S32 S2
TN32    L   KF  TN  S32 S2
HY64    L   KFC HY  S64 S1
NST64   L   KFC NS  S64 S1
TN64    L   KFC TN  S64 S1
HY65    L   KFC HY  S65 S2
NST65   L   KFC NS  S65 S2
TN65    L   KFC TN  S65 S2
HY19    L   NF  HY  S19 S1
NST19   L   NF  NS  S19 S1
TN19    L   NF  TN  S19 S1
HY24    L   NF  HY  S24 S2
NST24   L   NF  NS  S24 S2
TN24    L   NF  TN  S24 S2
HY05    S   ISO HY  S05 S1
NST5    S   ISO NS  S05 S1
TN05    S   ISO TN  S05 S1
HY12    S   ISO HY  S12 S2
NST12   S   ISO NS  S12 S2
TN12    S   ISO TN  S12 S2
HY31    S   KF  HY  S31 S1
NST31   S   KF  NS  S31 S1
TN31    S   KF  TN  S31 S1
HY34    S   KF  HY  S34 S2
NST34   S   KF  NS  S34 S2
TN34    S   KF  TN  S34 S2
HY62    S   KFC HY  S62 S1
NST62   S   KFC NS  S62 S1
TN62    S   KFC TN  S62 S1
HY63    S   KFC HY  S63 S2
NST63   S   KFC NS  S63 S2
TN63    S   KFC TN  S63 S2
HY04    S   NF  HY  S04 S1
NST4    S   NF  NS  S04 S1
TN04    S   NF  TN  S04 S1
HY20    S   NF  HY  S20 S2
NST20   S   NF  NS  S20 S2
TN20    S   NF  TN  S20 S2

My first attempt was building simple full (m1) and reduced (m2) models that gets directly at our question of interest but doesn't control for nested individuals.
m1 <- model.matrix(~ Region + Social * Injection, colData_filt)
m2 <- model.matrix(~ Region + Social + Injection, colData_filt)

We want to control for individual/batch effects, which is nested within both "Injection" and "Social" but not region, as we have three brain regions per individual. I followed the example in the DESeq2 manual for creating a term (ind.n) distiguishing individuals nested within groups, but now I'm not sure how to create the full and reduced model given that I have one more level than the example.
I've tried a really elaborate full model (m1) with the interaction term of interest (Injection:Social) removed for the reduced model (m2), but I'm not sure this is correct based on our design.
m1 <- model.matrix(~ Injection + Injection:ind.n + Injection:Social + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)
m2 <- model.matrix(~ Injection + Injection:ind.n + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)

I'm assuming this is wrong, but even if this was by some miracle the correct formulation, would there be a way to extract genes that explain the "Injection:Social" interaction term for separate brain regions?
As a work-around, I subsetted the data by region and ran three separate LRT analyses for each subset and compared the results. While this simplified the model to look like the first example above, I worry that we lose some power by ignoring the fact we have multiple brain region samples from single individuals across the dataset.
Any guidance is much appreciated. Thanks in advance

From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions.
So you set up the model like this:
m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..)

The last term should be Injection:Region and you can just use the waldTest (default) in DESeq2 for this term.
What does the terms do? ind.n*Region is the equivalent of ind.n + Region + ind.n:Region , and with this you effectively get an effect for every region in every individual.
Why don't we need the Injection:ind.n or Social:ind.n or Social:Region. These terms indicate the effect of Injection or Social can vary by individuals or regions. Most likely introducing too many parameters when you are interested in a common effect. Also you do not have the replicates or samples to distinguish this effect from region or other effects.
Since you provided an example we can run DESeq2 and you can see how the results look like:
mat = counts(makeExampleDESeqDataSet(n=1000,m=48))
dds = DESeqDataSetFromMatrix(mat,df,~ ind.n*Region + Injection + Social + Injection:Social)

dds = DESeq(dds)

resultsNames(dds)
 [1] "Intercept"            "ind.n_S2_vs_S1"       "Region_NS_vs_HY"
 [4] "Region_TN_vs_HY"      "Injection_S_vs_L"     "Social_KF_vs_ISO"
 [7] "Social_KFC_vs_ISO"    "Social_NF_vs_ISO"     "ind.nS2.RegionNS"
[10] "ind.nS2.RegionTN"     "InjectionS.SocialKF"  "InjectionS.SocialKFC"
[13] "InjectionS.SocialNF"

The terms you need are "InjectionS.SocialKF","InjectionS.SocialKFC", "InjectionS.SocialNF", and you can look at each of them:
head(results(dds,name="InjectionS.SocialNF"))
log2 fold change (MLE): InjectionS.SocialNF
Wald test p-value: InjectionS.SocialNF
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1  9.9811166787259   1.25304112986447 0.819806376919295   1.52845984752303
gene2 30.3449455820337 0.0329442893152027 0.705199688255367 0.0467162562092241
gene3 3.83223545055379    1.0281136369045  1.64095596190233  0.626533350543196
gene4  11.232305747171  0.595738624408923   0.8243883031544  0.722643227868976
gene5 6.70950627004097  0.756449993378065   1.0631622863378  0.711509430967263
gene6 26.1431134888287 -0.854784518963918 0.625714541243558  -1.36609342219393
                 pvalue              padj
              <numeric>         <numeric>
gene1 0.126398405431826 0.978671002658464
gene2 0.962739373909937 0.999897888026606
gene3 0.530965168963838 0.978671002658464
gene4 0.469899103018657 0.978671002658464
gene5 0.476768608734069 0.978671002658464
gene6 0.171909642630577 0.978671002658464

As mentioned, you can do a LRT if you want to test all the Injection:Social term interactions terms at one go, that is, the null hypothesis is that all of them are zero:
dds = nbinomLRT(dds,reduced=~ ind.n*Region + Injection + Social)
results(dds)

Usually the individual terms are more intuitively than testing all of them are zero, but you might have a special need for this.

Thank you, I will try this soon. Just to clarify, an LRT would still be appropriate for the larger question of finding the effect of the "Injection:Social" term and you suggest Wald for determining the region-specific impacts? Would it also be valid include a term like "Injection:Social:Region" in the suggested Wald model formulation in an attempt to explore the interactive effects of "Injection:Social" by region, or would that require more 3-way interaction terms with "ind.n"?

@jfaberha, thanks for providing the data, i have updated my answer, hopefully its clearer now. You can include "Injection:Social:Region" only if you include Region:Social. My point is this, do you think there is such an effect?

I would visualize this with a pca or correlation matrix to see where the effects are coming from before applying them onto the model. No need to have a super complicated model...

To answer your question regarding the "Injection:Social:Region" term, I will do some more research to see if this has an appreciable effect in the dataset. One question posed by the PI of this study is - what are the genes responsive to the "Injection:Social" interaction term in each of the different brain regions and are they different? I was guessing that maybe a 3-way interaction term could address this in the binomial Wald analysis, but I agree that adds excessive complexity to the model. That's said, I'm not sure how to get region-specific interactions without running separate analyses.

hmmm what is the main question? it sounds like everything under the sun. So it is indeed Injection:Social with different brain regions, you are better off running separate analysis for each region. It is easier for you to explain and honestly you don't gain power by analyzing all together

The main question is regarding the Injection:Social term across the dataset, with the brain region-specific interactions being less important follow-up. That said, I will move forward with LRT for the primary question and separate analyses for the region-specific analyses. Thank you for all the help!

sounds like a good plan!

Hi @jfaberha remember to upvote SW, LRTs are powerful statistics but never trivial to implement and that is a good answer.

in LRT can we use reference level which I use in pairwise comparison i have  a made a similar query Since i m not sure  if I use LRT which would be my reference in a reduced model


text:
What are the state-of-the-art cell-type RNA-Seq deconvolution methods?
I would like to find the proportion of each cell-type in bulk RNA-Seq transcriptomics data. I am looking for some guidance on the following:

What are the state-of-the-art methods?
What are their assumptions? (e.g. do they assume that cell-type expression signatures are available in advance?)

Can you describe your experimental design further? This is often solved with single cell RNA-seq

It's not my experiment - I am just trying to deconvolve bulk RNA-Seq data from a publicly available database (GTEx), e.g. given a whole blood sample, find out the proportions of each cell type (e.g. lymphocytes, macrophages,...). Any suggestion?

I see a lot of people using xCell. There should also be papers systematically compare many different methods, which might interest you.

Thanks for the pointer. Do you have any reference of a good up-to-date review paper? A simple description of the methods and requirements would be highly appreciated too

This one looks promising although it may not cover all the methods https://www.biorxiv.org/content/10.1101/2020.01.10.897116v1.full.pdf


text:
Do you know a program/script to count how many sequences have mismatches on each position of my primers?
I'm working with primer analysis and design, and so far I've performing it in Excel but: a) I don't like Excel, b) I do it manually, c) I don't like Excel, d) it gets tedious when working with thousands of sequences and tens of primers. My analysis consists of establishing how many sequences have mismatches on each position of my primers, and also count how many mismatches each sequence has with the corresponding primers (see attached image). This last analysis I've been able to (partially) solve it by using ncbi-blast+ tool. However, I will find very helpful a program/script where you input your sequences and primers fasta files and retrieve the number of mismatches at every primer position. Does anyone know something that could help me achieve all this? If it also has an output tab-separated file as the one attached that would be perfect for easy visualization of results. Suggestion/orientations are also grateful! I just want to know if there's something to help me out, otherwise I'll try writing a script to automatize my analysis. Thank you!

Its important to note this isnt a very clear question, for example how many loci are you examining, <10 or 100s? We tend to like knowing the organism and goal of the research. Anyway, at a guess I would just use NCBI's primer-blast program, here https://www.ncbi.nlm.nih.gov/tools/primer-blast/ . This appears to do everything on your list of requests

Michael, thanks for your comments! At the time, I'm working with >50000 sequences and 5 pair of primers/probes for one gene. Primer-blast does not support IUPAC nucleotide nomenclature, so that's a first limitation I've found. Also, I don't find the output information from primer-blast very helpful since it does not summarize the number of mismatches per primer position, and that's what I need in order to improve primer/probe design (I need a count of sequences rather than a list of them). Or perhaps I'm not setting correctly the output format parameters. Any further suggestions? Bests!

I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.
The hammer I use at the moment for this sort of stuff is LAST, which will do local alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.
I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are many other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:
Option 1 - Web-based LAST
In this, I use the Web version of LAST to generate alignments for the primers [as "reference"] mapped to the sequences [as "query"].

To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the Advanced Parameters section, set the a: gap existence cost to 1, and the b: gap extension cost to 1.
Click on Align!
Once the Alignment Results appear, click on the Alignments in tabular format link down the bottom, and save the file to your computer as last.tab.csv (the .csv at the end is the most important bit, to make sure that Excel recognises it as something it can open).
Open that file in Excel, set the delimiter to <tab>
The first column should be the alignment score, which in this case is the match score. This can be converted into a mismatch score by subtracting it from the primer sequence length (column 6).

Option 2 - Local LAST Installation
Installation
Install LAST and samtools.
Index Generation
The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the query, and the primer sequences will be the target. I use Martin Frith's lastdb mantra to do this:
$ lastdb -uNEAR -R01 primers.fa primers.fa

Mapping
The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to tab, and save the result to last.tab.csv:
$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

Aggregating data
Finally, it's a matter of aggregating the count information. I use grep -v to exclude lines that begin with a #, use awk to subtract field 1 from field 6, then sort the results and use uniq -c to count the mismatches:
$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2


text:
awk field separtor
I have a file like this :

and I want the output as:

I tried using awk :
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }'

and it prints the first id ENST.
But when I change it to:
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $5); print }'

it doesnt work.
What is wrong that I am doing here. Kindly help.

could you please provide a text-format toy example (or first few lines) of your input file? that would help to quickly see an application of your tests with awk.

done.KIndly see

thanks, although i meant more in an easily copy-pastable text format (not a screenshot) in a code snippet for example, same you used to write in your awk lines

Assuming you want to only show the fifth id from the first field, the following Perl script should work:
perl -pe 's/^(.*?)\t//; @F=split("\\|", $1); print(@F[4]."\t")' input.tsv > output.tsv

Explanation of the Perl code:

Loop over the lines in the input file (this is enabled by the command-line option -p)
Remove everything up to the first tab character from the input line, and save the removed part (excluding tab character) into the grouping variable $1 [identified by brackets].
Split variable $1 at the pipe character, |, saving the split components into the array @F
Print the fifth field [0,1,2,3,4] from the array @F [without a line break at the end]
Print the remainder of the input line (this is also enabled by the command-line option -p)

DDX11L1-202 is in the very first line of the excel document 5th character in the pipe.

You need the first column's 6th component when split by | and then the remaining columns of the file where \t is the delimiter. Use awk.
awk -F "\t" -vOFS="\t" 'NR==1 {print} NR>1{split($1,arr,"|"); print arr[6],$2,$3,$4,$5}' your_file

You could just do this in libreoffice calc using text to columns (data menu)
If you want to use awk (or maybe gawk) you could use capture groups; but rather that mess around with a fairly complex regex I'd probably do it in two stages like below
awk -F"|" '{print $5,$0}' OFS="\t" your_file|\
awk -F"\t" '{print $1,substr($0,index($0,$3))}'

The second awk prints everything except field 2

DDX11L1-202 1657|lncRNA| it retains the other columns with pipes.


text:
Question: Bismark methylation for non-uniquely mapping BS-seq reads
I am looking for the DNA methylation, specifically inside small RNAs; piRNA. I am mapping BS-seq reads on the genome with bismark. From the user guide of bismark, it maps reads uniquely on the genome, which mean it reports first read with best score and discard rest of the alignment and if a read mapping to multiple location with same score it will discard all such reads.
Typical command for obtaining uniquely mapped reads from bismark
~/Software/Bismark-0.22.3/bismark GenomePATH/Original_Genome/  --bowtie2 -p 10 -1 Methylation_reads/Sample1_R1.fq.gz -2 Methylation_reads/Sample1_R2.fq.gz --bam -o Sample1

Since, I am interested in piRNA only which may have multiple source in the genome (one piRNA coded from different location of genome), I wanted to take all mapped locations in consideration. Can you suggest what changes I have to make in bismark command. After reading bismark user guide I am confused because I could not able to find relevant parameter.
Also, please let me know if this approach is correct, if not what may be best approach for this.

I'm not certain about your piRNA approach, maybe others will be able to help with that question.
As for multi-mapped reads, from Bismark's documentation you can find the --ambig-bam flag:

For reads that have multiple alignments a random alignment is written out to a special file ending in .ambiguous.bam. The alignments are in Bowtie2 format and do not any contain Bismark specific entries such as the methylation call etc. These ambiguous BAM files are intended to be used as coverage estimators for variant callers.

It doesn't look like Bismark will return all of the multiple alignments, just a randomly selected one for a given read or read pair (this probably results from the way Bowtie2 handles multiple alignments within Bismark).

Thanks for the reply James. Does this ambig files can be used to call methylation just like normal bismark bam file using bismark_methylation_extractor? I read the manual and it seems that they cant be use for calling methylation.

As far as I can tell, it doesn't look Bismark adds the XM tags that would contain the methylation calls. So no, I don't think you can extract methylation calls from those reads, which is probably the safest option


text:
How do I color selective residues on the surface of a protein in pymol?
I am somewhat familiar with the pymol userface, but not very good with commands. I'm working with a pentameric subunit protein and I want to color specific residues on the surface of the protein, not just the cartoon. I can use select resi #-#, and color those with the (sele) tab on the right, but I do not know how to color the surface this way. When you select the sequence it will only color the cartoon. Any suggestions?
I'm using PyMOL 2.3.5 on Windows.

To clarify, you are seeing reds and blues on the surface? If so, these are the sidechain oxygens and nitrogens. So you need to colour the whole residue not only the element C atoms. color 0xFF7F50, resi 123 will colour the residue in the cartoon and surface as I omitted and element C. To revert color atomic.
In PyMOL an atom has a colour, which is then used by the representations. To have a differentially coloured representation, you need to make a copy of the object. create copied, obj_name and hide cartoon and show surface will do it.


text:
Modify how a script prints the output?
I'm trying to print out all measured distances between a ligand and 2 residues (using PyMol), below is the script that I'm using. I only want to measure the distances between all the ligand's O & ND2 of N312, as well as all the ligand's N & OD2 of D113.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "/r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("dst4.txt", 'w')

for a in stored.a:
  for b in stored.b:
    for c in stored.c:
      for d in stored.d:
        f.write( "distance from %s to %s is %s\n" % (a, c,
cmd.distance( "id %s" % a, "id %s" % c,cutoff=5,mode=2)))
        f.write( "distance from %s to %s is %s\n" % (b, d,
cmd.distance( "id %s" % b, "id %s" % d,cutoff=5,mode=2)))

f.close()

The output file currently looks like this
distance from 9 to 1334 is 4.589359760284424
distance from 8 to 4204 is 3.1509475708007812
distance from 9 to 1334 is 4.589359760284424
distance from 15 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 17 to 4204 is 0.0
distance from 9 to 1334 is 4.589359760284424
distance from 20 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 8 to 4204 is 3.1509475708007812
distance from 23 to 1334 is 0.0
distance from 15 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 17 to 4204 is 0.0
distance from 23 to 1334 is 0.0
distance from 20 to 4204 is 0.0

I want my output to look more like this. Thanks in advance!
Lig_N D113 4.59
Lig_O N312 3.15

Let's simplify...
cmd.distance(xx,xxx) accepts two selector strings and is fine with and state 1. And best of all if you select multiple atoms it gives you the smallest distance!
So all you need to do is:
# measure
d_O = cmd.distance('UNL and element O', 'resi 312 and name ND2')
d_N = cmd.distance('UNL and element N', 'resi 113 and name OD2')
# write
f=open("dst4.txt", 'w')
f.write(f'Lig_N D113 {d_N:.2f}\n')
f.write(f'Lig_O N312 {d_O:.2f}\n')
f.close()

If you are using python 2 (if so do note pymol scripting is better in 3), then the second bit is
f=open("dst4.txt", 'w')
f.write('Lig_N D113 %.2f\n' % (d_N,))
f.write(f'Lig_O N312 %.2f\n' % (d_O,))
f.close()

Unsolicited coding tips
A few unsolicited tips... Your original script is rather over-complicated and obfuscated and uncommented, so for better legibility why not do:

cmd.iterate_state is always harder to follow than the for atom in get_model('xxxx').atoms synthax
unless you are stuck with python 2, use f-strings not the %s sprintf style formatting
the /x/x/x/x selector synthax is never as readable as the longer
resi/resn/name selection algebra.

Code golf is always tempting, but it can result in future-you getting stuck as happened here...

You're right about the unnecessary complication in my script haha... I adapted the script a bit and it works :D Thank you!


text:
Determining significance of a variable in a glm model
How would one determine the significance of a variable in a glm model?
If I, for example, have a dataframe like the one seen below, how would I determine if the origin of the sample has a significant effect on the value? (this is the number of enzymes capable of degrading the substrate of that matter)
Substrate    variable value origin
cellulose       M09    8    free
mannan          M12    2    free
glycogen        M65    2    free
chitin          M87    4    free
cellulose       M90    2    isolate
manan           M78    1    isolate
glycogen        M21    4    isolate
chitin          M21    1    isolate

So far I have tried:
mcomp = glm.nb(value ~ origin, data = my_data)

summary(mcomp)
Deviance Residuals:
    Min       1Q   Median       3Q      Max
-0.9625  -0.9047  -0.9047   0.1212   3.5232

Coefficients:
              Estimate Std. Error z value Pr(>|z|)
(Intercept)   -0.01657    0.06571  -0.252  0.80097
originisolate -0.21911    0.08180  -2.679  0.00739 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(0.3418) family taken to be 1)

    Null deviance: 2053.5  on 2679  degrees of freedom
Residual deviance: 2046.3  on 2678  degrees of freedom
AIC: 6517.5

Number of Fisher Scoring iterations: 1

              Theta:  0.3418
          Std. Err.:  0.0186

 2 x log-likelihood:  -6511.4590

So free becomes the intercept and then isolate if significantly different from that. Does this mean Origin has a significant effect on the value?
Would the better approach be to do the following?:
mcomp = glm.nb(value ~ origin + Substrate, data = comb_data)
summary(aov(mcomp))
              Df Sum Sq Mean Sq F value Pr(>F)
origin         1     23   22.55   6.612 0.0102 *
Substrate     44   1445   32.84   9.631 <2e-16 ***
Residuals   2634   8981    3.41
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’

This shows me that origin and substrate have an effect on value if I understand correctly?

Second viewing of the question from what I can see -0.22 as a coefficient of origin is a strong negative association, so yeah it has a major impact. Its not how I would have done it, but that looks to be the result.

First viewing,
I'm going to throw my hat in here. We don't know what 'origin' is about, anyway just throw everything, i.e. each substrate and the origin into the same regression calculation. Check for a low-residual and preferably do a Q-Q plot, transform your data it this doesn't look good.
The key and the thing you are missing is your regression weights, without that I couldn't say very much. If the regression weight is near zero for 'origin' then it has zero impact. If the regression weight of 'origin' is positively greater than everything else ... I assume there are skewed distributions of 'substrates' between the 'origins'. If the regression weight of 'origin' is negative but still greater than all other regression weights then it is adversely affecting the 'value' you are seeking.
I don't know the experiment, the biological system or really the 'substrate' assays, so I can't comment any further.
The two issues I have are:

Doing an ANOVA on the output of a regression analysis doesn't make much sense to me. It is not something I would do, nor something in ML or GLM I've encountered.
Are you doing pairwise substrate/origin calculations? I presue not, but just in case this not how GLM works.

oh I was under the impression aov(summary(glm_model)) was a pretty common practice as it nicely summarises down to a variable level? Where do I find the regression weights?

When transforming the data do I have to transform both variables or can I transform just one of them until they fit the same distribution? and given that I am using a negative binomial glm does a uniform distribution between the variables matter as much?
I am unsure if there is a skewed distribution of variables between origins, that is what I am looking to find out really.

hi @Michael, very nice points here, yes the term anova is quite misleading sometimes. Normally people associate associate with the classical analysis of variance. In more "modern" uses, it can refer to the breakdown of deviance in a model, https://stats.stackexchange.com/questions/289226/what-are-the-differences-between-anovas-and-glms

It is like a log likelihood ratio test essentially, the "F" statistic instead of Chi-sq is used.

@StupidWolf So it is fine to aov(summary(glm_model)) to get a nice overview of the model? And by pairwise substrate/origin calculations does you mean doing the model on a subsatrate of the dataframe for each substrate @michael?

There is no better method, it's a matter of what you want to test or what is your question.
Using the anova() or aov(), test the terms collectively. For example, in your example with Substrate, the null hypothesis is that the coefficients are all zero, meaning cellulose =0, mannan =0 , ....
If the question is, "do the isolate samples have a higher value than origin samples?", then you can use your first model, where free is set as the reference and you test whether the effect of isolate is non-zero. Likewise you can do this for substrate and set of them as your reference. You can also do other pairwise comparisons using this model.
If the question is, "does origin have a significant effect on value, after controlling for substrate?", then you can use your second model.

Would I control for substrate by forming the model as such: value ~ origin + substrate ?

yes you do. you are giving every substrate a meaning  and on top of that estimating the effect of origin

And does the order in which I give the variable to the model, i.e. origin before substrate, indicate some sort of predefined importance to the model?

it doesn't define the importance, you can put it in any order. However if they are correlated, most of the variance will go to the first. I don't think this will apply in your case

see https://cran.r-project.org/doc/contrib/Faraway-PRA.pdf blocking design. A good ieda to read this

Thank you I will look at this :)


text:
kraken2 output to OTU table
I ran kraken2 on a 16s rRNA sample. I get two kinds of output, one like this
and the other one like this. So I have many files like this for each sample of my control and test group.
How do I generate OTU table out of these outputs for phyloseq usage? Do I parse the taxa ID column of any of these and generate a frequency table of each taxa in sample to generate OTU?
Any help or suggestion would be highly appreciated.

able to fix it ..issue resolved

Can you share how you fixed this please? Having the same learning curve!

Can you please add your fix / resolution as a separate answer (and accept it)? This is appearing in the question list as unanswered.

okay may be i will close the question..


text:
Calculate Entropy for DNA Multiple Sequence Alignment in R
I am pretty new to R. So I apologize for asking maybe a very basic question.
Let's say I have a fasta file with sequence below:
>sequence_1
ACCTGC--A
>sequence_2
ACC-GCTTA
>sequence_3
ACCTGCTTA

Is there a function or method to calculate entropy for each column in R?
I did some research and found a function that calculate entropy for Amino Acids but not for DNA nucleotides.
I would like to have the output be stored in a vector or in a list.
Thanks for the help in advance!

Hi what is the goal of your research? Entropy is a function of a given allele since you have variation, particularly indels this is not possible. With an indel position e.g. 4 now becomes position 3 (sequence_2) within an entropy calculation, so which column are you referirng to? I think you are confusing this with phylogenetic methods.

Hi Michael, I was able to get it to with some other people's help. But to answer your question, The MolecularEntropy function in R did exactly what wanted to. It calculated the entropy at each nucleotide position with aligned sequences.

okay, the blurb does say "Gaps are ignored on a site by site basis so the entropy values may have different number of observations among sites." It also needs sequences to be the same length ... I personally don't think it works with indels


text:
Relaxed and sequential Phylip format conversion
I routinely use AlignIO in BioPython for manipulating alignments and in context here moving from fasta format to phylip format. Alot of phylogeny packages accept a relaxed (sequence ID) & sequential (sequence) phylip format. This format is user-friendly, helps double-check for alignment errors and can be the essential input format.
A typical bit of code is,
'''
Converts a fasta file to phylip format
'''
from Bio import AlignIO

Class Convert (inheritClass):

    def fa2phy (self):
        fileOut, fileIn = self._inOut() # builds the path name and deletes any pre-exising outfile
        align = AlignIO.read(fileIn, "fasta") # reads the alignment in fasta
        AlignIO.write(align, fileOut, "phylip-relaxed") # writes in phylip

The output of "phylip-relaxed" is interleaved format
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAA

ACGTGACG .... ATCCCCAA
ACGTGACG .... ATCGGGAA
ACGTGACG .... ATCAAAAA
ACGTGACG .... ATTTTTAA

Output sought is relaxed-sequential
4 120
virusA_detailed_description ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_detailed_description ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_detailed_description ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_detailed_description ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

There is an option to allow a sequential format (above). HOWEVER, a "phylip-sequential" option output is as follows (it truncates the sequence ID to 10 characters,
4 120
virusA_det ACGTGACG .... ATCCCCAAACGTGACG .... ATCCCCAA
virusB_det ACGTGACG .... ATCGGGAAACGTGACG .... ATCGGGAA
virusC_det ACGTGACG .... ATCAAAAAACGTGACG .... ATCAAAAA
virusD_det ACGTGACG .... ATTTTTAAACGTGACG .... ATTTTTAA

No matter how I look through AlignIO options the above is the nearest format.  Could anyone point to a IO style solution?
I could write it, but that defeats the object (i.e. purpose not an instance of a method).

Note the "...." in the outputs above are not a typical "sequence identity" symbols, it just means that total length of a line in an interleaved format is 60 nucleotides (it might be 80).

I think the solution is to output using SeqIO, i.e. one sequence a time

I have nutted out a solution for you.
If this is the content of input.fa:
>virusA_detailed_description
ACGTGACGATCCCCAAACGTGACGATCCCCAA
>virusB_detailed_description
ACGTGACGATCGGGAAACGTGACGATCGGGAA
>virusC_detailed_description
ACGTGACGATCAAAAAACGTGACGATCAAAAA
>virusD_detailed_description
ACGTGACGATTTTTAAACGTGACGATTTTTAA

And this is the script:
#!/usr/bin/env python3

'''
Converts a fasta file to phylip format
'''

from Bio import AlignIO
from Bio.AlignIO.PhylipIO import PhylipWriter
import sys

infile = sys.argv[1]
outfile = sys.argv[2]
idwidth = int(sys.argv[3])
with open(infile, 'r') as input_handle:
    align = AlignIO.read(input_handle, "fasta")
    with open(outfile, 'w') as output_handle:
        PhylipWriter(output_handle).write_alignment(align, id_width=idwidth)

Then, save this script as phylpwrtr.py and do:
chmod +x phylpwrtr.py
./phylpwrtr.py input.fa output_longids.phy 30 # 30 is the width of the header, adjust the value to adjust header length

The contents of output_longids.phy looks like:
 4 32
virusA_detailed_description    ACGTGACGAT CCCCAAACGT GACGATCCCC AA
virusB_detailed_description    ACGTGACGAT CGGGAAACGT GACGATCGGG AA
virusC_detailed_description    ACGTGACGAT CAAAAAACGT GACGATCAAA AA
virusD_detailed_description    ACGTGACGAT TTTTAAACGT GACGATTTTT AA

Wow! Thank you  PhylipWriter is the solution ... huge thanks! :-)


text:
Is it possible to predict protein-ligand binding kinetics using machine learning?
I would like to work on a project that involves the prediction of protein-ligand binding kinetics. What might be the feature that is relevant for the prediction?

Does this answer your question? What is the current state-of-the-art in protein-ligand binding affinity prediction using machine learning?

Anyway, yes but you need access to a large database of experimental daa for protein-ligand affinity binding for training and validation. In addition, you'll need to vectorise the ligands ... which may not be at all easy. Personally, I would update your original question.

This is a comment to redirect to other answers, but go too long.

Docking is the in silico prediction of where a ligand binds in a protein.
There are several Q&As about docking here that are relevant, so please check out:

machine learning, docking and the factors that affect it: Generate ligands candidates based on protein shape —also worth a read the "learning from the ligand" paper paper, which discusses why some ligand only parameters tell a lot about binding.
correlation between docking scores and enzyme kinetics: Purpose of correlating docking score with catalytic efficiency(Kcat/km) —discusses a paper https://www.nature.com/articles/ncomms12965 wherein a linear regression is done between docking scores and kcat/KM


text:
R not being able to find the file despite being in correct directory
I have question regarding a project I have to do about a health student study. I cannot import the csv file via read.csv or read.csv2 despite having transformed all files to csv and being in the correct folder.
Here are the command lines:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")
getwd()
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

And here are the results in the console window:
setwd ("C:/Users/fek/Desktop/Analyse i-Share")

getwd()

[1] "C:/Users/fek/Desktop/Analyse i-Share"
TA<-read.csv2(file=Biobanque.csv, header=TRUE, sep=",", dec=",", fill=TRUE)

[1] Error in read.table(file = file, header = header, sep = sep, quote = quote,  :
objet 'Biobanque.csv' introuvable
To translate from french, it says it can't find the file!

The problem is just what it says.  By omitting the quotes around "Biobanque.csv", R thinks you are talking about an object whose name is Biobanque.csv. And R is confused because of course you have not created an object by that name.

This is correct - the error message states that R cannot find the object 'Biobanque.csv', not that it cannot find the file 'Biobanque.csv' (which would give the error "No such file or directory").


text:
Decontaminating RNA seqs for de novo transcriptome assembly and annotation of novel eukaryotes
I have raw paired-end RNA-seq reads for two novel eukaryotic species. Some background: the reads represent a copepod (arthropod) species each. The mRNA for each read set was obtained by extracting mRNA from multiple individuals. I am not exactly sure what sequencer was used (some Illumina short read machine as far as I can tell; I am not directly involved in this project).
I'd like to assemble these reads into transcriptomes de novo as no reference genomes are available as far as I know; my hope is to assemble with Trinity and then annotate with Dammit using its default databases (and check for completeness w/ BUSCO).
However, it looks like something went wrong during the sequencing runs (or during sample preparation) because the sequence GC content plots from FastQC look like this:

Here the curve with the tall peak is one sample, and the more flatter curve is the other. There are two closely related curves per sample because FastQC apparently seems to have recognized each read (forward and reverse) as a separate sample.
I first presumed that this is just rRNA contamination (I believe the mRNA was prepared via poly-A selection), and I filtered the reads against all of sortmerna's default databases (which I believe include sequences from SILVA and Rfam among others). Then I passed the sortmerna'd read through FastP to filter out adapters, paired read correction, and low complexity filtering (alongside whatever else FastP does by default). I had ran FastQC before FastP to check if the rRNA filtering did anything and it hadn't, so I also asked FastP to filter by length (--length_required 151) since there were small amounts of reads with lengths shorter than this value which I thought could be the source of the additional peaks in the data.
However, even after sortmerna and FastP together, nothing has really changed:

(There are only two traces/curves here instead of four because the paired reads were apparently recognized properly by FastQC now.)
My plan was to run Khmer to clean the reads a little bit more and then assemble them. I did do a quick assembly with Trinity (with the raw reads) and the BUSCO scores were pretty bad (only ~ 25% of the arthropod BUSCO genes present in the assemblies). I was hoping I'd be able to get Dammit to annotate whatever was left after the cleaning, despite the poor BUSCO prognosis.
Now I'm not entirely sure what to do if I can't clean the reads to begin with, since preventing nonsensical contigs would be absolutely important for a decent annotation.
So this leads me to my questions:

Are these data even salvageable?
Is there anything else I can do to clean the data?
Is this even contamination to begin with that I'm seeing here? (I am presuming some kind of contamination because both samples share that peak at ca. 18% GC content, among other shared features.) Would mapping all the reads against each other and discarding the shared sequences help? If yes, is there a tool I can use for this?
Any other thoughts in this regard?

I searched through the posts here and on Biostars as well as I could but I wasn't really able to find anything helpful in this regard (hence my own post).
I'm very new to all of this, so any inputs, tips, advice, and/or guidance would be most appreciated.

Please can you provide the other metrics given by fastQC? eg the quality per base? I didn't use the latest version by i am suprised by the y-axis in %. In addition, how did you run Trinity? on a merged set of reads or did you perform 1 assembly for 1 sample?

I am sorry but the links to fastqcfiles seems broken. Can you just run a fastp (without filtering according to length, or set a lower size threshold, eg 36 nt) followed by Trinity? You may also give a try to Bridger or BinPacker which are faster than Trinity without quality loss.

also, can you please tell me how you got the data. An academic platform? a private company?

according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them.
Try to run fastp with the following options:
fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

and run FastQc again to check to efficacity of the trimming.

sorry I can't see them on the repo... but can you provide the json or html output of fastp too? thanks!

ok thanks. Indeed, no strong improvement. I think you may provide custom fasta sequences to fastp using the --adapter_fasta argument. You may add sequences from Clontech Universal Primer and Clontech SMART CDS kits because contamination most likely come from that. Besides, do you have an idea about the number of genes in the genomes you are interested in? I saw you had less than 10M of reads and i wonder if it is sufficient for a complete de novo transcriptome assembly.

yes you're definitively! Concerning the overrepresented sequences, I never faced this problem for now. But the % is quite high (more than 3%) so i guess it deserved to be tested. As an alternative, you may try to map your RNA-seq reads on the related genome species (use HiSAT2+Stringtie pipeline). It might help to clean your reads!

for sure, do not hesitate to contact me again. Hope you will find a solution for your data anyway.


text:
Difference between isoforms and paralogs in transcriptoms?
I've assembled RNAseq data into a transcriptoms using Trinity. There's a option to keep only the longest isoforms for each transcripts and it lead me to wonder how it deals with duplicated genes (-->transcripts). I haven't found an anwser in the original paper. Althought I reckon paralogs and isoforms are two distinct features I don't see a way to make a distinction when you only have trasncripts.

It depends. If you are performing Illumina based RNAseq, which is likely, it is definitely tricky, You should be able to do this in an ideal world because,

a paralogue will contain SNPs, i.e. there will be nucleotide variation between paralogues, but the mRNA can of an identical length
isoform will also have identical sequence but the length of the mRNA will vary depending on what initiation codon it used

I sort of agree however because,

there will be 'grey-zones' where paralogues are different sizes and VCF/ SNP calling is not perfect and ...
assembly of paralogues can be difficult, risking collapsing a paralogue into a single mRNA, but it is not impossible depending on the level of SNP variation
assembly of isoforms would be really hard

I personally wouldn't know how point 3 is solvable, but I guess it is if your isoform is relatively small, i.e. well within the read-length of Illumina output.

If a reference genome assembly exists, then reference-guided assembly with Trinity will keep paralogs, because it treats each genomic region as a different assembly set.
If a reference genome doesn't exist, paralogs will probably result in the construction of non-representative "isoforms" containing variants present in multiple genes. Inchworm does a greedy assembly of reads, so if a chimeric transcript can be constructed from the reads (which is likely with paralogous genes), then it might end up being fed to the subsequent steps. Whether that is actually emitted at the end by Butterfly will mostly depend on how the paths through the putative assemblies are traversed, which probably depends on the tolerance of Bowtie to errors.
If paralogous transcripts are sufficiently similar that the differences are within the error tolerance for assembly, then I would expect that they'd be collapsed into a single transcript.

Thanks @gringer its worth doing assemblies both ways, IMO it would work well if your 100% on all paralogue and isoform possibilities. In organisms with massive multigene families (where you don't know the full repertoire of sequences) using de novo assemblies is doable to prevent the collapse of paralogues depending on the distance between mutations and is particulalry true now Illumina now has pretty long read lengths. In the days of the 150 -200 bp reads then yes collapose would happen.  Clearly minority paralogues can be lost regardless.


text:
How can I obtain a list of all NCBI gene ID's along with their full name, symbol, and also known as symbols?
I would like to do two things:
(1) Obtain a list of all NCBI (Entrez) gene ID's.
(2) Obtain the "Official Symbol", "Official Full Name", and "Also known as" fields from the NCBI gene pages for each of those genes.
This information, since it is displayed on the NCBI website, must be located in some databases, but I have not had luck finding a tool that will allow me to extract this information.

All of these data are in the gene_info.gz file on the NCBI FTP site: https://ftp.ncbi.nlm.nih.gov/gene/DATA/
Specifically, you may want to look at the gene_info.gz file. It has the following fields (shown for an example gene below):
                              #tax_id [  1]: 9606
                               GeneID [  2]: 5768
                               Symbol [  3]: QSOX1
                             LocusTag [  4]: -
                             Synonyms [  5]: Q6|QSCN6
                              dbXrefs [  6]: MIM:603120|HGNC:HGNC:9756|Ensembl:ENSG00000116260
                           chromosome [  7]: 1
                         map_location [  8]: 1q25.2
                          description [  9]: quiescin sulfhydryl oxidase 1
                         type_of_gene [ 10]: protein-coding
   Symbol_from_nomenclature_authority [ 11]: QSOX1
Full_name_from_nomenclature_authority [ 12]: quiescin sulfhydryl oxidase 1
                  Nomenclature_status [ 13]: O
                   Other_designations [ 14]: sulfhydryl oxidase 1|quiescin Q6 sulfhydryl oxidase 1|testis tissue sperm-binding protein Li 62n|thiol oxidase 1
                    Modification_date [ 15]: 20200604
                         Feature_type [ 16]: -

If you are starting  with a set of NCBI GeneIDs and want this sort of information for just those, you can use Entrez Direct for that. If you are starting with RefSeq transcript or protein accessions, be sure to check out the other files in the FTP path mentioned above.

I had tried using this file, but a bug in my code made me think it didn't have what I wanted.  Your answer made me look into it again, and I found the bug.  It's exactly what I needed.

Go to Ensembl (https://www.ensembl.org/index.html), click BioMart on the menu at the top of the page. Select the dataset you want to use, then go to Attributes, select the feature of each gene you want to get. Then click "Results" button to make sure that the table looks like what you want. Then output the results to a file.

This is almost exactly what I'm looking for.  The only thing it is missing is as comprehensive a "Also known as" list.  For example, NCBI lists the PPP4C gene has having a total of 7 symbols, while this data set has 2.


text:
Multiple input & output files?
I'm trying to print out the measured distances for each .pdbqt file, and then write a new output file for each .pdbqt, using the original name of the .pdbqt file as the new .txt file name. Below is the script.
I'm using python3 & I'm trying to do this via Linux CLI.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("tt.txt", 'w')

for a in stored.a:
  for c in stored.c:
      d1=cmd.distance("id %s" % a, "id %s" % c,cutoff=5)
      f.write("%s D113 %s\n" % (a,d1))
      print(d1)

for b in stored.b:
  for d in stored.d:
      d2 = cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
      f.write("%s N312 %s\n" % (b,d2))
      print(d2)

f.close()

My folder looks something like this
.
├── dst6.txt
├── sp
│   ├── r55_01.pdbqt
│   ├── r55_02.pdbqt
│   ├── r55_03.pdbqt
│   ├── r55_04.pdbqt
│   ├── r55_05.pdbqt
│   ├── r55_06.pdbqt
│   ├── r55_07.pdbqt
│   ├── r55_08.pdbqt
│   ├── r55_09.pdbqt
│   ├── r55_10.pdbqt
│   ├── r56_01.pdbqt
│   ├── r56_02.pdbqt
│   ├── r56_03.pdbqt
│   ├── r56_04.pdbqt
│   ├── r56_05.pdbqt
│   ├── r56_06.pdbqt
│   ├── r56_07.pdbqt
│   ├── r56_08.pdbqt
│   ├── r56_09.pdbqt
│   └── r56_10.pdbqt
└── t6.py

THANKS IN ADVANCE!!

You made the hardest work, you were not far to achieve what you wanted ;)
This script.py should work after you made it executable (in a terminal type: chmod u+x script.py), just type ./script.py in your terminal in your working directory (so here you should be in your "sp" directory). Note that I got inspired by: https://www.biostars.org/p/42474/ to know how to launch pymol within the script and also that I would personally not use pymol for doing this kind of analysis/data extraction from a pdb/pdbqt file but rather use a dedicated analysis package (there are plenty like MDanalysis). I personally use pymol (pymol script) for visualization purpose.
#!/usr/bin/env python

import __main__
__main__.pymol_argv = ['pymol', '-qc']  # Quiet and no GUI

from pymol import stored
import pymol, os

pymol.finish_launching()

#put what you've done into a function
def myfunc(file):
    spath = os.path.abspath(file)
    sname = spath.split('/')[-1].split('.')[0]

    pymol.cmd.load(spath, sname)

    stored.a, stored.b, stored.c, stored.d = [], [], [], []
    pymol.cmd.iterate_state(1, "%s///UNL`1/N*" % sname, "stored.a.append(ID)")
    pymol.cmd.iterate_state(1, "%s///UNL`1/O*" % sname, "stored.b.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASP`113/OD2" % sname, "stored.c.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASN`312/ND2" % sname, "stored.d.append(ID)")

    f = open('%s.txt' % sname, 'w')

    for a in stored.a:
        for c in stored.c:
            d1 = pymol.cmd.distance("id %s" % a, "id %s" % c, cutoff=5)
            f.write("%s D113 %s\n" % (a, d1))
            print(d1)

    for b in stored.b:
        for d in stored.d:
            d2 = pymol.cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
            f.write("%s N312 %s\n" % (b, d2))
            print(d2)

    f.close()

mypdbqt = os.popen('ls *.pdbqt').read().split()

for i in mypdbqt:
    myfunc(i)

pymol.cmd.quit()


text:
Can SBOL encode DNA Assembly linkers?
Where in an SBOL file could you specify the linkers that would be placed as prefixes / suffixes between parts for assembly? Lets say you have a final construct of your design and want to assemble the component parts - could you encode what linkers you used between each component, ie as an annotation? Would this be different between SBOL versions?

Every version of SBOL allows you to specify prefixes/suffixes as annotations on a sequence.

In SBOL 3, the construct would be a Component of type DNA (https://identifiers.org/SBO:0000251), and the linker a SequenceFeature on that Component. The role of the sequence feature might depend on the particular assembly you are describing, and whether you're talking before or after assembly (in those cases where the process changes the sequence), but a good one in many cases would be Nucleotide Binding Site (https://identifiers.org/SO:0001655)
In SBOL 2, the representation is the same, except the construct would be a ComponentDefinition and the linker a SequenceAnnotation.
In SBOL 1, the representation is still the same, except the construct would be a DnaComponent and the linker an SBOL 1 SequenceAnnotation (which has different fields than in SBOL 2).

If the assembly method does not change the sequence, then all three can represent this pretty much equally well.  If you want to represent a sequence that is modified by assembly or the plasmid carriers from which parts are extracted, then SBOL 3 is the best because it cleanly represents deletions during assembly via the sourceLocation field.


text:
Seurat analysis only with matrix with no barcode and features
I have downloaded a dataset from GEO and only the matrix (mx) files are available. How can I perform the analysis without the barcode and the gene features?

what do you mean by gene features? Does the matrix have genes as rownames? Can you share the GEO number? you need gene names or metadata ...

Seurat has different functions that can read different types of data. Please share a few rows and columns of your data and we can try to point you to the right direction (function).


text:
Do I need to fill in missing side chains and loops before protein-ligand docking for VLS?
I am going to screen a set of small molecules against a CYP450 2C9 protein target. The structure I am using exactly is 5A5I. Do I need to fill in missing side chains and loops before performing VLS? I assume while preparing the protein I should leave water molecules in the structure as there is a bond between $Fe^{2+}$ and water that is involved in the binding of a ligand.

In general
Most if not all docking applications will tolerate a jump.
Furthermore, you really don't want to dock anything to the flexible regions even if you are using a program that does properly protein flexibility.
Furthermore, if your protein adopts a closed conformation upon native substrate binding (not really the case here), you should start with that for better docking results.
Cytochrome P450
So most Cytochrome P450s are strange though and in many those loops play a large part in their catalysis. I am on a paper under review (bioarxiv) that deals exactly with this in a different cytochrome (P450-BM3, the biotech one). Namely even if the active site is fairly large for an active site in the core of an enzyme, the ligand (testosterone) did not really rotate around and instead the specificity of oxidation was controlled by how the ligand got in. So those crystallographically unresolved loops act as gatekeepers for what can bind. In terms of VLS, one aims for blocking off the cofactor, i.e. the heme. So you can do docking without worrying about them in a first instance, but there is a possibility that a lot of active drugs instead interact with the loops.
Certainly, dynamic undocking and other MD methods may be required a few iterations in, to see why some candidates failed and why others worked.
An interesting unofficial further possible spin is the fact that solvent-exposed sequential residues Asp-Pro can hydrolyse under acidic conditions, in a mechanism similar to isoaspartate formation, and several cytochromes have these on the entrance loops and have weird kinetics and smeary SDS-PAGE gels.
Lastly, you are totally correct in your observation. The water binding the heme in the crystal is very special. It is taking the place of oxygen. That could be modelled as is, changed to a pre-encounter O2, or event changing the heme to a heme bound to a radical oxygen —note that some forcefields (CHARMM, Rosetta Ref2015) do not like this weird valency iron.


text:
how to extract a different substring for each row of a dataframe in R?
I have a dataframe with putative alignments of thousands of probe sequences. I would like to pull out the adjacent nucleotide from the genome for each alignment. My dataframe includes seq.name for the target chromosome and I've added an adjacent_pos variable with the location on the target chromosome that I want. Using read.fasta from the phylotools package gives me a genome in the format of a dataframe with seq.name and seq.text for each chromosome and unplaced scaffold.
Working on a test data set with a single chromosome, I successfully used something like this:
adj <- cbind(alignments$adjacent_pos, alignments$adjacent_pos)
alignments$adjacent_base <- str_sub(genome$seq.text, adj)

Now back to the whole genome, I can't wrap my head around how to direct str_sub to the correct seq.name for each seq.position. I'm sure there's an elegant solution, but I haven't been able to find it!

please can you provide a small reproducible example for your alignments and genome objects? the merge() function from the data.table package will surely help you.

Say your file looks like this:
cat test.fa
>chr01
ACCGGTT
>chr02
ATGCATGC

Something like what you have done, to read in the fasta and extract the sequence. Key part is creating a GenomicRanges object:
library(GenomicRanges)
library(BSgenome)
library(Biostrings)

adj <- data.frame(chr=c("chr01","chr02"),start=c(2,3),end=c(2,3))

    chr start end
1 chr01     2   2
2 chr02     3   3

adj = makeGRangesFromDataFrame(adj)

fasta <- readDNAStringSet("test.fa")
SEQ = getSeq(fasta,makeGRangesFromDataFrame(adj))

as.character(SEQ)
[1] "C" "G"

or if the fasta file is huge, you can use Rsamtools, indexing the fasta file first:
library(Rsamtools)
fasta = FaFile("test.fa")
indexFa('test.fa')

as.character(getSeq(FaFile("test.fa"),adj))
chr01 chr02
  "C"   "G"

but it's not just one location from each chromosome. I have about 160k alignments, each needing a different location pulled out from a 3G genome that has 21 chromosomes plus over 7500 unplaced scaffolds. Still perhaps the GenomicRanges object is the way to go?

yes.. the code above works for as many locations as you have. You just need to construct the genomic ranges object like I did above, and I suggest using the Rsamtools solution if your file is huge

adj should be GenomicRanges object...if it was clear from the solution above

thank you I will try

brilliant, worked like a charm! thank you

you're welcome :) glad it worked for you


text:
Error after trimming illumina adapters
I am removing illumina adapters of the NGS data with a loop. My NGS data is storage in /data/HTS_seq/.
I used this function:
 for infile in /data/HTS_seq/*_read1.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

 for infile in /data/HTS_seq/*_read2.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

and It gives to me this error:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'
  Traceback (most recent call last):
  File "/usr/local/bin/cutadapt", line 10, in <module>
cutadapt.main()
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 652, in main
options.discard_untrimmed)
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 327, in trimmed_and_untrimmed_files
trimmed = untrimmed = xopen(output_path, 'w')
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/xopen.py", line 178, in xopen
return open(filename, mode)
      IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which is my error?
Thank you in advance.

First of all, you aren't actually setting the $outfile variable. Your loop has this basic format:
for var in *; do newVar="something" command $newVar; done

For example:
$ for var in 1 2 3 4; do newVar="new:$var" echo "$var: $newVar"; done
1:
2:
3:
4:

As you can see, newVar is never actually being set. That's because setting a variable is its own operation. So you need to do that, end it with a control character like ;, or &, or a newline, and then run your cutadapt:
$ for var in 1 2 3 4; do newVar="new:$var"; echo "$var: $newVar"; done
1: new:1
2: new:2
3: new:3
4: new:4

Or, with your data:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

But that would give a weird file name, like foo_S1_L001_read1.fastq_trim.fastq. You probably want names like foo_S1_L001_read1_trim.fastq instead, so you would want:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Now, the specific error you show makes me think that isn't the command you ran. Could you perhaps have run this one instead by mistake (note the / instead of \ before the _trim.fastq)?
for infile in /data/HTS_seq/*_read1.fastq; do
 outfile=$infile/_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Your error message is:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which means that cutadapt tried to open /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/ as a directory and then look for a file named _trim.fastq inside it. The only way I can think that this would happen is if you ran the command I gave above and not the one you show in your question.

It works! Thank you Terdon

There is a bug here
/data/HTS_seq/*_read1.fastq

This will go to your root directory / and there is definiately no folder called data in the root. To prove it type ls / and that is where the command is pointing, then compare this to ls ./.
To solve it go to the directory HTS_seq and type pwd then copy and paste the directory structure. I suspect ./data/HTS_seq will solve it
 for infile in ./data/HTS_seq; ...

This will look for the folder data in your current working directory rather than the root directory.

Okay you've been assigned /data as root ... no biggie.
PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq

Why is / in the file name? What happens if you ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq?
and ls /data/HTS_seq/ Finally could you please ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq\/_trim.fastq?
The \/ is important. What I think has happened is a / is in the file name and that will cause problems because it will think that _trim.fastq is a directory when its part of a file name, and would explain the error.
Then do the same thing using head instead of ls see what happens.
If the last ls and head works then rename the files by typing rename in the command line. If rename works it will be easy to remove / from all your files using the power of the rename function (its a Perl function BTW).

Hi Michael, thanks you for your answer. I typed pwd in HTS_seq folder and I got this direction: /data/HTS_seq                                                                                Anyways I typed the code with the . :                                                                                for infile in ./data/HTS_seq/*_read1.fastq; do outfile=$infile_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done and I got the same error. Which is my mistake? Thanks you

Hi @AdriánP.L. please try the above, I should have read the error more carefully first time around

You can never have a / in a file name. The / alone with \0 are the only two characters that can never be in file names.


text:
How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?
I am using edgeR to perform differential expression (DE) analysis on a set of RNA-seq data samples (2 controls; 8 treatments). To correct for batch effects, I am using RUVSeq.
I am able to get a list of DE genes without normalization:
x <- as.factor(rep(c("Ctl","Inf"),c(2,8)))
set <- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, "DE_genelist.txt", sep=""))

Then immediately after creating the "top" object, I use RUVg to normalize:
# [...]
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg <- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, "DE_RUVg_genelist.txt", sep=""))

And I get the error:
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("SeqExpressionSet", package = "EDASeq")’ to a data.frame

I am not sure how to print the list of normalized results like I can with the unnormalized data. Ideally, I would get a file with the same format as the edgeR output (as a .csv or .txt file):
"logFC" "logCPM" "LR" "PValue" "FDR"
"COBLL1" -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
"UBE2D1" -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
"NEK7" -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
"SMC6" -2.300 5.674738981329 61.8130019860261 3.7767230643666e-15 3.4974443325016e-12

How can I get a list of genes as an output after normalization with RUVSeq?

You do the normalization before running your edgeR. The purpose of RUVg is to remove "Remove Unwanted Variation Using Control Genes". In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.
Using the example dataset in vignette:
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set <- betweenLaneNormalization(set, which="upper")

set1 <- RUVg(set, spikes, k=1)

You can look at it, it's an expression set with counts etc, not results:
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples
  element names: counts, normalizedCounts, offset
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

You run edgeR now on the results of RUVg:
design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)

This is great! Definitely answers the question. [Although the very last line topTags(lrt) will only provide a list of 10 genes, so I'll use topTags(lrt, n=nrow(set))$table for a full list].

Follow-up question for you then: What is the utility of performing the "upper-quartile" (UQ) normalization before using RUVg and then running edgeR? My thought was, "if I am using RUVg to correct for batch effects, why should I also use UQ?" Perhaps it is important?

This is what I know about UQ normalization: scaling with or without UQ before RUVg will usually lead to similar results. It is only important to include the edgeR offset -- calcNormFactors(y, method="upperquartile") -- if you decide to scale prior to RUV use. If no UQ scaling is applied, the edgeR offset should be set to 0 (zero).

The betweenLaneNormalization() in limma adjust for signal intensity difference (in arrays) or in edgeR, you try to normalize between samples such that most genes are unchanged. RUVg adjust for unseen variation on top of this, for example batch effects. So you don't want to introduce this as a kind of variation for RUVg model to pick up. Hence you do the normalization before. And as you noted, if you use edgeR,  RUV also ask you to provide the offset , which is the normalization factors.

Hope this answers your question. You can try and see for yourself what is the offset, and what RUVg does on top of it

Wonderful explanation! Thanks for your help @StupidWolf!!

I have not used this package but from your code it seems that ruvg is not a table. Instead, it is an R object, which means that you cannot use write.table. I think the results you want is stored in the object. All R objects contain "slots" of data, which can be accessed by @. If I were you, I would type ruvg@ and should be able to see which data slots are contained in the object.

This certainly gives me more direction than I had before, so I appreciate that. Just looking at all the data stored by the ruvg object (Large SeqExpressionSet (1.1 Mb)) tells me that it may take some time to find what I am looking for in there. I'll report back if I find the results I want. Thanks for your suggestion!


text:
How do I convert an Excel file (.xlsx) to an .smi file format?
How do I convert an Excel file(.xlsx) to .smi file format?
The original file I have is an excel file containing over 27,000 ligands. I am going to using the Padel Descriptor software to calculate the descriptors for the ligands however, Padel Descriptor just supports .smi files. Screenshot of the excel file added below for more clarity.

So a smi file is simply a text file formatted SMILES tab name newline.
Therefore what you need to do is

save for safety
make a new sheet
copy the SMILES column and the name column
save as tab separated file
in Finder/Nautilus/File Explorer rename the file extension to .smi


text:
How do I get repeating regions from fasta file or directly from internet
If I go to ncbi site, apply customize view -> Customize -> All features -> Update view and then ctrl + f : Alu I can get all alu repeaing regions coordinates from this chromosome.
It looks like:
 repeat_region   complement(37490..37800)
                 /note="AluY"
                 /rpt_family="SINE/Alu"
                 /rpt_type=DISPERSED

How could I download all this coordinates and names of repeating regions from the internet to .txt file? Or some tool

Hi @WhoAmI, an efetch command would be the best option. You can Blast them out ... but the genomic position might not follow (not sure).

@Michael could you write usage example?

Have a go yourself, thats the best way to get help on site. I don't work on human genomics so its not something I'd do ... but fishing out bits of Genbank annotation is certainly an efetch thing.

Click on the 'GenBank' dropdown menu at the top left corner of this page to see a graphical view of the sequence. You should see a 'Repeat region' track in the Sequence Viewer. Click on the 'Download' button in the top right corner of the Sequence Viewer, choose 'Download Track Data' to launch a dialog box where you can enter a range, select the 'Repeat region' track and download data in BED format. Unfortunately, the names for the features are not great.
If you want the names to be the same ones in the flat file, you can download the feature table and parse it. To download the feature table, click on the 'Send To' link at the top right corner of this page, choose 'Complete Record', 'File' and 'Feature Table' as 'Format'. After you download the feature table, you can parse it using awk as shown below to create a 4-column table with the data.
cat seq.ft \
  | awk 'BEGIN{FS="\t";OFS="\t"}{
    if(NF==3) {print s,e,f,q ; s=$1; e=$2; f=$3; q=""}
    else if(NF==2) {print s,e,f,q ; s=$1; e=$2; q=""}
    else if(NF==5) {q=q"|"$4"="$5} }' \
  > feat_table.tsv


text:
Selecting part of an extracted ligand
I have the PDB 6LU7.

I extracted the ligand using extract lig, org  and then the protein using extract prot, poly.  The ligand has three parts:  02J, PJE, and 010.

I want to select the 02J and 010 IN A SCRIPT and delete them.

I know how to do it by simply clicking, but want to incluse this in a script.

Hi @user8338 and welcome. I think what you are asking is PyMol scripting, I suspect Modeller will do something similar. Do you use PyMol/Modeller? If so could you please update and clarify your post accordingly. Thanks!

I know that structure way too well —SARS-COV-2 MPro.
In PyMOL command line
 remove resn 02J or 010

or
 cmd.remove('resn 02J or 010')

In python
 import pymol2
 with pymol2.PyMOL() as pymol:
      pymol.cmd.fetch('6LU7')
      pymol.cmd.remove('resn 02J or 010')
      pymol.cmd.save('6LU7.mod.pdb')
      #or pdbblock = pymol.cmd.get_pdbstr()

From the sounds of it you'd like to make the peptidomimetic a peptide. If you wanted to do the horrid thing of making PJE, the capped side chain glutamine ketamide, into glutamine, the command to rename atoms is a multitude of cmd.alter('resn PJE and name <atomname>', 'name="<newname>") followed by a cmd.alter('resn PJE', 'resn="GLN"') and lastly cmd.sort() to commit the changes. Not pretty, but is all in PyMOL.

Wow, impressive response - thanks

And, related, for the structure 6M2N, there are four 3WL ligands, one in each of the four protein chains.  I'd like to remove 3WL in chains B, C, and D, retain chain A ligand, also scripted (not select and remove).  Thanks!

The selection algebra is listed in https://pymolwiki.org/index.php/Selection_Algebra
But basically you want to select based on meeting two criteria so it's an and. Namely, cmd.remove('3WL and not chain A')

Typo: cmd.remove('resn 3WL and not chain A')


text:
gene set enrichment in single-cell data
Gene set or pathway analysis in scRNA-seq data has its own challenges. Since the data is sparse and many genes will be missing for any given cell, some people simply add up all the counts across the gene set. The sum is treated as the expression value of the gene set. I still see this simple approach used in recent papers. Seurat AddModuleScore (and used by CellCycleScoring) is basically the same method, except it additionally subtracts the signal from a control gene set, but that is just shifting the score down.
My concern is that the gene set score can be driven almost entirely by a single gene. If the max value of one gene in a gene set is 100x greater than for others, which certainly happens, then those other ones are not really contributing. What is even the point of using a whole gene set then? Yet, this approach is widely used, so it must be working well enough. What am I missing?


text:
Programmatically retrieve Accession List from NCBI
I would like to programmatically retrieve the file "Accession List" from an NCBI page.
I'm running firefox, so I'm clicking the file and clicking Ctrl + Shift + E to try and figure out what the file endpoint is. But I can't seem to figure it out (I haven't done this much, I just remember doing it once for a previous project).
When I do this, I see two links popping up on my screen, as shown below:

I see that the first one has a field (I'm talking about the panel on the right) called "acc", which varies from one data dump to the next, so I'm presuming this is the link I need to use to retrieve the file (though I'm not sure), and just change the "acc" field depending on the data I want. Assuming this is right, I'm not sure how to go from here.
I've included the text on the right below for clarity.
GET
    https://www.google-analytics.com/collect?v=1&_v=j83&a=1402353855&t=event&_s=4&dl=https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA491191&o=acc_s%3Aa&ul=en-us&de=UTF-8&dt=Run Selector :: NCBI&sd=24-bit&sr=1280x720&vp=1232x601&je=0&ec=Link Click&ea=Metadata  &el=NONE&_u=SACAAEAB~&jid=&gjid=&cid=1350895533.1590937485&tid=UA-68648221-3&_gid=391407234.1594002343&cd1=run_selector&cd2=sra&cd3=&cd5=NO&cd6=run_selector:sra:&cd8=&cd9=&cd10=&cd11=&cd12=&cd13=&cd14=&cd15=&cd16=&cg1=run_selector&cg2=sra&cg3=&cg4=run_selector:sra:&z=508671370

Which accessions are you trying to get to? In the page you refer to, there is an option to download the list of accessions - one SRR run accession per line - by clicking the 'Accessions' button or the entire table by clicking the 'Metadata' button as highlighted in the screenshot below:

You can get the same data from the command line using the Entrez Direct. The command for that would be something like:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo

to download the table in CSV format.


text:
GSEA enrichr with 10x genomics differential_expression ranks
I am attempting to use GSEA enrichr with 10x genomics differential_expression rankings. Reading what people seem to be doing with GSEA, there seems to be a pre-/post-singlecell gene expression versioning of the tools associated, so I want to make sure what I am doing makes sense.
Explanation:
For each gene in columns 1 and 2, there are 3 columns in each of the clusters that the 10X Genomics cellranger count software produces, e.g.:
Cluster 3 Mean Counts   Cluster 3 Log2 fold change   Cluster 3 Adjusted p value
0.2392551160479023      -6.591931364165717           0.029177998608334674
0.1435530696287414      -5.608385347948343           0.0003833420249831818
1.107409394278862       -5.441829814465051           3.8059567638748296e-09
0.1948220230675776      -4.9155704967404             0.008739337761029243
0.20849374398460058     -4.65702722396748            1.0980978990257375e-11
0.36571853453036496     -4.6498340596466585          1.1896434946397147e-17
0.1811503021505546      -4.445327449243366           1.3416233207430882e-12
0.4409129995739914      -4.433552335619653           1.2400396349891802e-16
1.7226368355448967      -4.240367394494225           2.2072280062348652e-15
0.10937376733618391     -3.902215199856403           2.0315030904651526e-12
3.2948847410025404      -3.7580328996321235          0.028858857004902847
0.3588826740718535      -3.738237134795285           8.439629283346646e-13
[...]

In trying to give a tissue / cell type identity to each of these clusters, I've used the following method to feed a subset of genes to EnrichR:
1- Sort the Log Fold change column numerically so we get low negative values.
2- Filter out for the Adjusted p-value to have a value smaller than 0.05.
3- Filter out for the Means Counts value to be bigger than 0.01.
4- Take the top n genes from the list, and feed it to GSEA EnrichR to compare against Mouse_Gene_Atlas or ARCHS4_Tissues.
Given that the fold change of the clusters is relative to each other, does it make more sense to rank the top n genes numerically or in inverse numerical value? E.g. selecting the first 100 genes in each cluster that have the most negative Log2 fold change values or the most positive Log2 fold change values or both (the top 100 genes with highest absolute Log2 fold change)?
Looking at previous answers, I get that one shouldn't rank by p-value: e.g. see @gringer reply
https://bioinformatics.stackexchange.com/a/13465/180

tagging @gringer who authored an answer to a related question.

You should be using the most positive fold changes, not the most negative. The positive ones are the ones that are most highly expressed in your cluster of interest. You probably want to use less than 100 genes. Depending on the experiment, you may not have that many good markers for a particular cluster.
The other question you mentioned refers to DESeq2 specifically, so it is not applicable here.

What's a good number smaller than 100? 50? Thx

Depends on your dataset. It's usually very difficult to have 100 genes that identify the cluster well. You can pick the last gene and see how well it actually marks that cluster.

would it be adequate to trim down the list of genes by making sure they are not present in other clusters? E.g. find the X-most genes in Cluster 1 where none or only, say, 10-20% of them are in Custers 2 to N. Is that a reasonable approach? Thanks in advance.

What do you mean by "not present"? Do you mean detectable in only one cluster? If they are only detectable in one, then they should already be your top cluster markers by fold change or p-value.


text:
Combining several different affymetrix expression microarray experiments for further MDS analysis
Assume that I have three different Affymetrix microarray experiments, 2x U133 (from 2015 and 2017) and 1xU133 Plus 2.0 (2019). Each array contains several different samples, some of them are common for 2 or 3 arrays. I'd like to combine these thee results to perform further MDS analysis using whole sets.
I'm aware about some batch effect that I'd like to avoid.
When I do MDS using separately normalized data from three experiments, I can see that expression pattern of sample present in two of 3 matrices is the same but the normalized signal values are raised by a certain value. It visibly affects further MDS analysis.
My questions are:

Is there some bioconductor package or strategy to do such analysis avoiding bias and batch effect?
Can I just normalize all 3 experiments together using RMA? And then MDS analysis will be unbiased.

Thanks in advance!