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

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

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

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

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

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

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


text:
DMPfold: number of models produced using "10 500" feature in "run_dmpfold.sh"?
I ran this command:
run_dmpfold.sh sequence.fasta sequence.21c sequence.map ./folder 10 500

From my previous discussion, I understand that I should run with default "run_dmpfold.sh" parameters i.e. 3 50; however, since I have already run the 10 500 I have some questions.
I see that after generating "sequence_1.pdb to sequence_500.pdb" files, these files are no more in the folder and, now I see "ensemble_1.pdb".
Will there be more ensemble files like "ensemble_500.pdb" and "sequence_n.pdb" files? I mean when can I expect "final_1.pdb' file? Also, I see files like "CLUSTERS_001.pdb".

The sequence_n.pdb files are temporary and are combined into an ensemble.i.pdb file after each iteration i. That file has many structures in. After the run is complete there will be a final_1.pdb file. Don't worry about the cluster files. If you are interested in what they do, consult the run_dmpfold.sh script and see what it is doing line by line.
As I have said multiple times elsewhere, you should run the script with default parameters and let it run to completion. Then use the final_1.pdb and any other final_*.pdb files as the output. If you want more structures, use the ensemble.i.pdb files.

source: Stack Exchange

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

sort: votes


text:
How I deal with this kind of gene expression comparision
I have RNA-seq from two sequencing batches; Lab technician says that he has run the RNA expression quantification two times in bathes 1 and 2 for example tumor 1 in batch 1 and tumor 1 in batch 2 , normal 2 in batch1 and normal 2 in batch 2. PCA shows however shows two experimental runs are too close to each other and can be merged
[![enter image description here][1]][1]
Now my question is can I compare only a tumour sample with its own matched normal sample in one patient to get differentially expressed genes for this patient by using experimental runs (technical replicates )? Design would be Tumour1_batch1 Tumour1_batch2 versus Normal1_batch1 Normal1_bacth2
> head(mycols)
              Condition
Normal1_batch1  N
Normal1_batch2  N
Tumour1_batch1  T
Tumour1_batch2  T

Yes, you could (although I think you would need three samples in each side to be able to calculate the sd and var). However, how generalizable it would be ? What is your biological question behind this comparison? One thing is to pair the samples and control for batch effect and another is to just compare 4 samples.

My main question here is: what are the differentially expressed genes between tumour and normal adjacent sample for this specific patient where I am going to use two experimental runs as my replications for having statistical power  . Thank you

Using technical replicates tells you how good your technician can use a pipette, it gives you no information about the biology behind it. In statistics we want to sample (take samples from a population) to get info about the whole population. For example, if you want to have the average length from a population, you sample different people to get a good average. In your case you would measure one person twice, it tells you nothing about the population but it tells how good your measuring skills are.

If you want to compare individual patients, why are there less normal samples compared to tumor samples?

Because I don't have real replication rather two experimental runs

No. If the only difference between Tumour_batch1 and Tumour_batch2 was the library prep/sequencing run, then all this will tell you is if the particular sample of RNA taken from the tumour is different from the particular sample of RNA taken from normal. It almost certainly is.
Indeed if you did technical replicates of two RNA samples, both taken from normal, you would almost certainly find differences.
What you want to find is the differences in RNA level between tumour in general and normal in general. To do that you need to make an estimate of the mean and variance of RNA levels in the tumour and in the normal. To do that you need multiple samples from each.
If you need differences within a single patient, it would be more honest just to compute the fold changes without any implication of statistical certainty. You would then need to think about alternate schemes for measuring the confidence in inferences drawn from the data when you consider the actual biological questions.

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


text:
Parsing specific features from Genbank by label?
I'm trying to parse a Genbank file to find a specific feature. I can pull it out if I know the feature type (e.g. repeat_region). E.g. if I'm looking for this feature:
 repeat_region   5623..5756
                 /label=5' ITR
                 /note="5' ITR"

I know that I can find it using:
for feature in reference.features:
    if feature.type == "repeat_region":
        print(feature.location)

But I don't trust that it will always be a repeat_region. Instead, I'd like to look for it by label (5' ITR). I can seem to find a way to parse that from the feature object. Any suggestions?

The /label=5'ITR is called a qualifier. You can look for those qualifiers:
for feature in reference.features:
    for k, v in feature.qualifiers.items():
        if k == "label" and v == ["5' ITR"]:  # or "5' ITR" in v
            print(feature.location)


text:
1000 genomes missing SNPs on chr X in GRCh38
I am trying to get a reference together for imputation using the 1000 genomes data set. The GRCh38 version of the data set was made available earlier this year, but seems to be missing a significant number of SNPs on the X chromosome (~96%).
Is there any resource that can explain why? I have tried searching on their home page with no luck.
The data sets in question are GRCh38, GRCh37.

Please do not cross post to BioStars, Stack Exchange and 1000 Genomes helpdesk as it can create duplicate work:
https://www.biostars.org/p/399254/

As I told you on 1000 genomes helpdesk:
Because of the ploidy settings we used for our variant calling, we only called variants in the the pseudoautosomal region on the X chromosome. We are working on settings to correctly analyse X variants in both males and females and plan to analyse X properly soon.

I think these are biallelic snps:

The set is restricted to biallelic SNVs and INDELs.

If you check the liftover version, the chromosome X is complete...


text:
Biopython Genbank writer not splitting long lines
I am parsing a csv file of annotated sequences and using Biopython to generate Genbank files for each. I want to add annotations of the sequence features.
My output file shows features listed without the correct line breaks. Other software is then unable to parse the names of the features
misc_feature    58..60  /note="Oxidation"
misc_feature    100..102  /note="Oxidation"
misc_feature    241..243  /note="Oxidation"

It should be formatted like this
misc_feature    58..60
                /note="Oxidation"
misc_feature    100..102
                /note="Oxidation"
misc_feature    241..243
                /note="Oxidation"

My section of code adds a space between the location and the /note="" but the Biopython writer doesn't add line breaks.
for i in range(len(start_list)):
    feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
    feature.key = "misc_feature"
    liab_type = "{}".format(liability).split(":")[1]
#need to add 1 to start list since csv counting is off
    feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i]) + "  " + '/note=\"{}\"'.format(liab_type)
    container.features.append(feature)
    i +=1

The Biopython source is toggling the wrapping of the feature qualifiers on no_space_key and Bio.GenBank._BaseGenBankConsumer.remove_space_keys
How do I enable wrapping of text to format output correctly?
    def __str__(self):
    """Return feature as a GenBank format string."""
    output = Record.INTERNAL_FEATURE_FORMAT % self.key
    output += _wrapped_genbank(self.location, Record.GB_FEATURE_INDENT,
                               split_char=',')
    for qualifier in self.qualifiers:
        output += " " * Record.GB_FEATURE_INDENT

        # determine whether we can wrap on spaces
        space_wrap = 1
        for no_space_key in \
                Bio.GenBank._BaseGenBankConsumer.remove_space_keys:
            if no_space_key in qualifier.key:
                space_wrap = 0

        output += _wrapped_genbank(qualifier.key + qualifier.value,
                                   Record.GB_FEATURE_INDENT, space_wrap)
    return output

Glad you found the solution!, The way to signal that the question is solved is posting your own question and selecting it as the one that resulted more helpful. This will unify the style for all questions and answers and make it easier for people to find the solution to your question.

I read the Biopython Genbank record source code carefully and realized I wasn’t using the Qualifier class correctly.
After reading the source code again, I realized I was using the data structure incorrectly.  I need to use the Qualifier class to add a key and feature.  This is formatted correctly and solves the problem.
feature = Feature()
# Strict Genbank format calls non-standard things misc_feature
feature.key = "misc_feature"
#need to add 1 to start list since csv counting is off
feature.location = "{}..{}".format(int(start_list[i]) +1, end_list[i])
liability_qualifier = Qualifier()
liability_qualifier.key = '/note="'
liability_qualifier.value = '\"{}\"'.format(liability).split(":")[1]
feature.qualifiers.append(liability_qualifier)

Thanks for posting back and explaining the answer to the rest of us


text:
Can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?
I am trying to perform de novo structure prediction of a 1135 long protein sequence; can DMPfold and I-Tasser work well for de novo structure prediction of 1135 amino acid long sequence?

For DMPfold you might like to look at Figure 5B of the paper, which implies that performance starts to drop off after ~600 residues. Only proteins up to 500 residues were used in the training set, for example.
A protein of that length is almost certainly multi-domain. I would consider doing domain prediction; we have a tool DomPred, but many others are available. Then run the sequences of the domains individually and recombine them later.
I don't know about I-Tasser's ability to do this, perhaps see if they mention it in their papers.

Generally, but not always the I-Tasser algorithm is not in the same league as a covariant mutation-based analysis (psicov) and don't think many people would disagree with this. I don't do I-Tasser but I'm aware it is the algorithm of choice for rough and ready predictions.
Psicov based models are essentially de novo predictions, based on amino acid mutation patterns, which I assume are refined against known homologous structure in DMFold (inaccuracies can creep through). I suspect the starting point of I-Tasser is a known homologue and it's associated known structure.
Having said that you need alot of sequence data to make psicov covariance work (>8000 aligned sequences) and (strictly in my opinion) it will also depend on the speed of mutation. Higher eukaryote  gene alignments could require much higher numbers of taxa in the alignment. The covariance value will give you the strength of the signal and values >0.7 should be used. I personally consider>0.5 to 0.7 borderline, which is fair enough because those are the generic thresholds in statistics. I don't know if you can access the underlying covariance values in DMFold, or moreover set the threshold.
So I-Tasser could outperform DMFold in given situations, but the underlying output of covariance based approaches should make it clear how well the approach has performed. Obviously if it performs well I'd go with that.
I would follow@jgreener 's advice on the follow up program.


text:
Changing the axis limits of ggplot objects
By EnhancedVolcano R package I have this plot with this code
    > head(res1)
         baseMean log2FoldChange     lfcSE        stat     pvalue      padj
SAMD11   237.7975     0.03584334 1.3244737  0.02706233 0.97841002 0.9971691
NOC2L   4751.0381    -0.62090622 0.6500115 -0.95522347 0.33946470 0.9239759
KLHL17   735.0849    -0.95126291 0.8417444 -1.13010892 0.25843033 0.8926461
PLEKHN1 1117.5598    -2.11321359 1.0411108 -2.02976816 0.04238011 0.6350143
PERM1    409.0061    -2.32444927 1.1255813 -2.06511000 0.03891259 0.6285916
HES4    1088.7549     0.28395931 0.9091259  0.31234320 0.75477971 0.9866881
>
p= EnhancedVolcano(res1,
    lab = rownames(res1),
    x = 'log2FoldChange',
    y = 'pvalue',
    xlim = c(-5, 8))

I have stored my plot in an object
> class(p)
[1] "gg"     "ggplot"

I want to limit the axis length in this plot in the other word shrining my plot more but I don't know how to do that. For example y axis be from 0 to 10 instead of 20
Any help please

Could you clarify what do you man by "limit the axis length"?  Do you want it for instance to reach just from -4 to 4 or do you want a wider plot? What kind of searches did you make? Also if you could include a reproducible example, see reprex, it would be easier for people to understand and work with your problem

Axes for ggplots can be controlled via xlim() and ylim() functions. You can try:
your_ggplot_object +
  ylim(-5,20) # assuming that you want the gene labels to be readabe, not clear from your post

p1 +

xlim() +

ylim(-5,20)

Error in UseMethod("limits") :
no applicable method for 'limits' applied to an object of class "NULL"

I have realized that the function EnhancedVolcano has a "ylim" argument and its lower limit is 0, please try "ylim = c(-5, max(-log10(toptable[,y]), na.rm=TRUE) + 5)" within EnhancedVolcano.

Sorry what is the y here? when I am placing pvalue, says that Error in [.data.frame(res_responders, , pvalue) :
object 'pvalue' not found

Please check the documentation of the function EnhancedVolcano() and play a little bit different arguments. For you case ssetting the ylim argument "ylim = (-5,20)" might help.


text:
How many models should be generated by rosetta for denovo protein structure prediction?
How many models should be generated by rosetta for de novo protein structure prediction of a 382 long sequence? In the tutorial, it says 50,000 models but practically it takes a lot of time (weeks) to generate this number. For instance, I am currently running rosetta on a 382 long sequence and in a week it has generated 1500 models; I am using 16 cpus with Openmpi.

Unfortunately you do in general increase your chance of success by going to the 50k range of models or higher. Projects like Rosetta@home exist because even large research groups can't access as much compute power as they'd like. In this case the long time arises from the combinatorial explosion of possible fragment combinations that need to be explored.
Success is rather protein dependent though - sometimes the right structure falls out after only a few models, other times your fragment library does not cover the protein adequately and you could run Rosetta forever without getting a good structure. Have a look at your 1,500 models and see how similar they look.
382 residues is rather long for Rosetta. Another factor to consider is whether you are using predicted contacts as part of the folding process (e.g. in the PconsFold protocol), which will greatly improve the chances of getting a good model and reduce the number of models you need. In fact the PconsFold group switched to CONFOLD for PconsFold2, I think mainly because it was much faster, so you are not alone in finding Rosetta runs long.

I did look at the generated 2100 decoys (rosetta ab initio run is still in progress) and plotted the RMSD vs ENERGY plot and it looks like this: https://ibb.co/LgQPz0L . Please let me know of your comments. I think there is still a long way to convergence? And, all structures are similar? Thanks

The RMSD is zero for all models, indicating something is wrong. Do you have the native structure for this protein? The RMSD on an RMSD v. energy plot should be between the the model and the native structure.

There is no native structure as this is an ab initio prediction. Could you please let me know about your comments?  Command line is this: "mpiexec -np 16 /main/source/bin/AbinitioRelax.mpi.linuxgccrelease -database /main/database/ @options -mpi_tracer_to_file log > log1.log &". "Options" file looks like this: https://ibb.co/V9VrQHV

I will answer where you have asked this as a question: https://bioinformatics.stackexchange.com/questions/10453/rosetta-ab-initio-run-gives-0-00-rmsd-for-generated-2000-models-run-still-in-p


text:
module differential connectivity (MDC) analysis
In this paper I am citing the lines below:

Intra-modular connectivity of resilient and susceptible network modules was compared to that of corresponding genes in control mice to identify gain, loss, or no change.Pie charts (lower panel) summarize module differential connectivity (MDC) analysis. Proportionally, very few resilient modules had significant MDC compared to more than half of all susceptible modules, which predominantly showed gain of connectivity
A module differential connectivity (MDC) metric quantified differences in coexpression networks to identify phenotype-specific modules
(Zhang et al., 2013)

Need to clarify whether this "module differential connectivity (MDC)" analysis is a part of WGCNA or its an independent analysis. Because if I do a search for module differential connectivity(MDC)analysis the result I get is this `DGCA: Differential Gene Correlation Analysis'. But no where in the paper they have mentioned "DGCA" library in the citation.
It would be helpful if i get some suggestion regarding this analysis.
This is the paper:
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3677161/

Answer from @devon-ryan, converted from comments:
They described how they did it in that paper, there's just not a tool for it. It's not part of anything. They describe the algorithm, maybe DGCA implemented it, maybe not, but they don't mention using it. WGCNA was obviously used to define the modules, though.


text:
Putting labels of different sizes on one PyMOl Object
I'm new to PyMol (and StackExchange!) and working on my first project. I have the structure of a protein as an object, called PolyA-M, and the idea is that the residues of it are shown as spheres of differing sizes based on a calcualted conservation value. I want to label each sphere with it's amino acid, but to have the label size corresponding to sphere size. Here is a sample of the code that might be used to label one amino acid.
alter ( resid 138 ), resn = "Q135"
alter ( name CB and resid 138 ), vdw = vdw * 0.8 * 0.679934640522875
set_color col138, [ 1.0, 0.181, 0.181 ]
show spheres, name CB and resi 138
show sticks, ( name CA or name CB ) and resid 138
set label_size, 14, PolyA-M
label PolyA-M and name CB and resid 138, "Q"
color col138, name CB and resi 138

In this case, residue 138 is medium sized sphere with a size 14 label Q that fits well. However, if later on in the code another residue is labelled as so:
alter ( resid 198 ), resn = "A208"
alter ( name CB and resid 198 ), vdw = vdw * 0.8 * 0.359803921568627
set_color col198, [ 1.0, 0.601, 0.601 ]
show spheres, name CB and resi 198
show sticks, ( name CA or name CB ) and resid 198
set label_size,  8, PolyA-M
label PolyA-M and name CB and resid 198, "A"
color col198, name CB and resi 198

This will label residue 198 with a small A to fit the small sphere. However, this will also change the earlier label to size 8 as well, making it too small for the sphere.
Is there any way for me to prevent this, such that each residue keeps its own label size?
I know that one way is to create multiple identical objects and keep all labels of the same size restricted to the same object (For example, PolyA-M_10 contains all residues with label size 10, polyA-M_14 all residues with label size 14) but I was wondeirng if there is a more efficient way?
Thank you for any help! This is my first Stack Exchange question so any feedback would be appreciated :)

I ended up fixing this by just putting this into the boilerplate code:
copy PolyA-MC_4, PolyA-MC
copy PolyA-MC_6, PolyA-MC
copy PolyA-MC_8, PolyA-MC
copy PolyA-MC_10, PolyA-MC
copy PolyA-MC_12, PolyA-MC
copy PolyA-MC_14, PolyA-MC
copy PolyA-MC_16, PolyA-MC
copy PolyA-MC_18, PolyA-MC
copy PolyA-MC_20, PolyA-MC

This copied my object such that I had one per size of label text, and then each label was placed on the object. As all labels on a given object are the same size, it worked! e.g. for one residue the labelling is now:
set label_size,  4, PolyA-MC_4

Effectively, it's the same as what I was already doing, but instead of having ten copy filies in the directory only one is needed. Much easier for sending to people, and more efficient!

'label_size' is a object-state-level setting, which means that you will have to rely on the 'create' command to create a new object for every different label size.
It can certainly be automated, and for this I'd recommend getting into Pyhton Scripting for PyMOL

Ah, I see what you mean! The way I've been doing it with multiple objects is to just create multiple copies in the directory (So PolyA-M_1 to PolyA-M_10) such that it reads those in too, and then writes labels of different sizes onto them. But if there isn't a command which is able to change label size at a residue level, then this definitely sounds more efficient than having loads of different files. It would also be good for me to learn more scripting in general. Thank you!


text:
Create a Heatmap that uses modules of genes rather than individual genes in Seurat
Is it possible to create a Heatmap in Seurat that takes in and displays modules of genes as its features rather than individual genes?  So if I have a module of genes associated with a trait or phenotype I can compare the expression of that module across clusters against the expression of other modules of genes.
Currently the way that I am trying to do it is similar to the code below (In this example I am only using one module to test if it works):
Module_28 <- readLines(con = "Path_to_file/module28.txt")
Module_28 <- list(Module_28)
Combined <- AddModuleScore(Combined, features = Module_28, ctrl = 5, name = 'Module_28')
names(x = Combined[[]])
avg <- AverageExpression(Combined, features = NULL, add.ident = NULL, return.seurat = TRUE, verbose = TRUE)
DoHeatmap(avg, features = c("Module_281"), group.by = "ident", size = 3, angle = 0, draw.lines = FALSE) + theme(axis.text.x = element_text(size = 7, angle = 45), axis.line = element_line(color = "#ffffff")) + scale_fill_gradient2(low = '#1000ff', mid = "#ffffff", hight = '#aa0101', space = "Lab", na.value = "#ffffff", midpoint = 0, guide = "colourbar", aesthetics = "fill") + coord_flip(expand = TRUE, clip = "on")

Often When I try to run the DoHeatmap function I get the following error:
Error in DoHeatmap(avg, features = c("Module_281"), group.by = "ident",  :
  No requested features found in the scale.data slot for the RNA assay.

What's wrong with your current approach?

I get an error that I'll add to my post above so you can see it.

Did you scale your data at some point?

@haci yes I did scale the data, with the ScaleData() function.


text:
How to convert VCF format file into Fasta(.fa) with the Reference data
I need a Software or Python Program for Converting VCF(.vcf) to  FASTA(.fa) format with the help of reference (.fa)

Not quite a dupe since it asks for an R or python approach, but very relevant: Converting a VCF into a FASTA given a reference with Python, R.

https://bioinformatics.stackexchange.com/questions/2825/converting-a-vcf-into-a-fasta-given-a-reference-with-python-r

Also vcftools..

I used FastaAlternateReferenceMaker previously. You will need to download GATK first. Next you prepare your reference genome (reference.fasta) and your vcf file (input.vcf). You will call the GATK function like this:
 java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -V input.vcf \

output.fasta will contain the new fasta with snps inserted at sites specific by the vcf file

I got this error when I used gatk on a Linux plaftorm:
A USER ERROR has occurred: Input files reference and features have incompatible contigs: No overlapping contigs found.

your chromosome names are different

how do I resolve this please?
Additionally, the vcf files are from somewhere external and I actually subset it to get a specific sample. my reference is https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz


text:
How to load and split my two individual datasets when integrating datasets in Seurat?
I want to do comparative analyses of my two individual single cell RNA-Seq datasets in Seurat. And I load and list my data as follows:
Control.data <- Read10X("~/Documents/...  Control_filtered_feature_bc_matrix")

Control <- CreateSeuratObject(counts  = Control.data, project = "Control_summary", min.cells = 3, min.features = 200)

Mut.data <- Read10X("~/Documents/... /Mut_filtered_feature_bc_matrix")

Mut <- CreateSeuratObject(counts  = Mut.data, project = "Mut_summary", min.cells = 3, min.features = 200)

Object<- list(Control,Mut)

names(Object) <- colnames(Object)

Then the protocol https://satijalab.org/seurat/v3.1/immune_alignment.html said, I need to use SplitObject function, but what do I split by?
Object_list <- SplitObject(Object, split.by = "?")

Can you please edit your post so as to include the link for the protocol that you have mentioned? Moreover, you have formatted your code as a "quote", please also change this to "code".

Seurat is expecting individual datasets to be normalized separately prior to data integration. In that respect lists of objects corresponding to different datasets are handy to manipulate each object/dataset individually.
In your case you do not need to do any splitting, you have an object corresponding to each of your datasets.
If you would have started with a single count matrix of all datasets then you would need to split your data based on batches/experimental condition/...


text:
A program "PSICOV" "stopped" from a few days on Ubuntu 18.04?
I am running dmpfold, "seq2maps.sh" and from a few days the program is at the "PSICOV" stage:
psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

"ps aux" gives me this:
daniel     26693  490  2.5 2546588 2514320 ?     Rl   سبت20 25116:58 psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

Any idea on what could be the matter? Yes, I am runinng "rosetta" too on my 32 CPUs, so all my CPUs are covered to max 100%; howver, RAM is 93Gb and only 32 GB has been used. Thanks

In summary just be patient. psicov doesn't  parrallelise so a tenfold increase in CPU will not hurry it along.
Also psicov is not RAM intensive given the machine you are describing, so you haven't hit a RAM bottleneck.
I don't know the ps output format you are using I'm use to seeing %CPU and %MEM from ps -aux on Susie and OSX. You should also try,
top -c -u username

or for OSX top -u
You should see the psicov happily chuntering away consuming ~100% of a single process. There is a UNIX/Linux guru specialist on BioSE who could almost certainly provide better diagnostics.


text:
Rosetta ab initio run gives 0.00 RMSD for generated 2000 models? (run still in progress)
I am running rosetta ab initio prediction using the command line:
"mpiexec -np 16 /main/source/bin/AbinitioRelax.mpi.linuxgccrelease -database /main/database/ @options -mpi_tracer_to_file log > log1.log &"

"Options" file looks like this: https://ibb.co/V9VrQHV
RMSD vs Energy plot is:
https://ibb.co/1mjzwDj
There is no native structure as this is an ab initio prediction. Could you please let me know about your comments on the RMSD being 0.00 for all 2000 models generated till now? Is this normal and will produce results after more runs?

RMSD is a measure of similarity between two structures. The relevant RMSD in protein structure prediction is between a model and the native structure. If you don't have a native structure you cannot have an RMSD in this sense. Ignore the fact it is 0.0, there is possibly some error that means you are calculating the RMSD between a model and itself.
And that's the difficulty of protein structure prediction... you don't know whether you are right or not! The structures you have generated are your results, and you'll need to find ways to filter or use them that are not reliant on RMSD. You could use estimation of model accuracy methods or clustering, for example.

thanks, when following the tutorial:https://www.rosettacommons.org/demos/latest/tutorials/denovo_structure_prediction/Denovo_structure_prediction; it says " However, if no related structure is known, you would use e.g. the lowest energy model (and rescore the top models with e.g. -native best_E.pdb)." Can you let me know how can I "rescore" the models? Thanks


text:
split 1 column input into 5 column bed file
In f1 below I am trying to split $1 based on each line and create a bed file. If the line is a snp then the : is spilt and the text is $1 the last digit is $2 -1 and $3 and the letter to the left of the > is $4 and the letter to the right is $5. (lines 1 and 2)
If the line has dup then the : is spilt and the text is $1 the last digit is $2 and $3 and the $4 is - and the letter to the right is $5. (line 4)
If the line has del with nothing else then the : is spilt and the text is $1 the last digit is $2 and $3 and the letter to the left of the > is $4 and $5 is -. (line 5)
I am not sure how to format for line 3 (the complex indel)... maybe split $1 on the : and the text is $1 the last digit is $2 and $3 and the letters to the left of the > is $4 and the letters to the right is $5. Thank you :)
f1
chr7:140453145A>T
chr7:140453136A>T
chr7:140453135_140453136delCAinsTT
chr20:31022287dupA
chr19:13054614delG

desired
chr7    140453144   140453145   A   T
chr7    140453135   140453136   A   T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

awk
awk 'BEGIN{OFS="\t"} {sub(/[^0-9]+$/, "", $1)}' f1

awk 'BEGIN { FS = "[ -]"; OFS="\t" } NF==3 { print $1, $2 - 1, $2, $3 } NF==4 { print $1, $2, $3, $4 }' f1

Do you have to do this in awk for some reason? Sure it's possible, but there are enough cases that'd it'd make more sense to use python that's easier to write longer scripts in.

python is fine I just don't know it as well.... only the basics but I am interested to see what it would look like as I can learn from it. Thank you :).

The following Python script seems to do the job.
#!/usr/bin/env python3

from __future__ import print_function
import re
import sys

for line in sys.stdin:
    if '_' in line:
        match = re.search(r'(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)', line)
        print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
    elif 'dup' in line or 'del' in line:
        match = re.search(r'(\S+):(\d+)(del|dup)([ACGT]+)', line)
        if match.group(3) == 'del':
            print(match.group(1), match.group(2), match.group(2), match.group(4), '-', sep='\t')
        else:
            print(match.group(1), match.group(2), match.group(2), '-', match.group(4), sep='\t')
    else:
        match = re.search(r'(\S+):(\d+)([ACGT]+)>([ACGT]+)', line)
        coord = int(match.group(2))
        print(match.group(1), coord - 1, coord, match.group(3), match.group(4), sep='\t')

Invoked like so on the command line.
[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

Thank you very much: print(match.group(1), match.group(2), match.group(3), match.group(5), match.group(4), sep='\t')
^
SyntaxError: invalid syntax
~$ python --version
Python 2.7.9  (not sure if this helps, sorry new to python. :)

It looks like you're running Python 2, which unfortunately is the default version on many operating systems even though it's at its end-of-life. I'll update the script so that it's compatible with Python 2.

"print" was a statement in Python 2 but is a function in Python 3 hence the syntax error.

The Python script in my first response is heavy on regex matching, which is pretty clunky in Python. I like Python much better than Perl overall, but a throwaway script like this will be clearer and more concise in Perl.
#!/usr/bin/env perl
use strict;

while(<STDIN>) {
    if (m/(\S+):(\d+)_(\d+)del([ACGT]+)ins([ACGT]+)/) {
        print("$1\t$2\t$3\t$5\t$4\n");
    }
    elsif (m/(\S+):(\d+)(del|dup)([ACGT]+)/) {
        if ($3 eq "del") {
            print("$1\t$2\t$2\t$4\t-\n");
        }
        else {
            print("$1\t$2\t$2\t-\t$4\n");
        }
    }
    elsif (m/(\S+):(\d+)([ACGT]+)>([ACGT]+)/) {
        my $coord2 = $2;
        my $coord1 = $2 - 1;
        print("$1\t$coord1\t$coord2\t$3\t$4\n")
    }
}

It's been several years since I wrote Perl on a regular basis, so there are probably ways to make the script even more clear and concise. (Although I have a love/hate relationship with the default/hidden variables, which can make a script more concise but also harder for a newcomer to understand.)
Invoked like so on the command line.
[standage@lappy ~] $ ./transform < f1
chr7 140453144 140453145 A T
chr7 140453135 140453136 A T
chr7    140453135   140453136   TT  CA
chr20   31022287    31022287    -   A
chr19   13054614    13054614    G   -

Quick Perl solution (following on from Daniel's answer saying that this would be "clearer and more concise in Perl"). I was a bit confused with your BED output because you've treated the start/end locations of a single base deletion and insertion different to a SNP. In my example here, the start and end locations for those cases are the same to distinguish from the two-base change for TT->CA [this can be changed]:
cat 10454.txt | \
perl -ne 'if(/^(.*?):([0-9]+)(_([0-9]+))?(del|dup)?([ACGT]+)((ins|>)([ACGT]+))?/){
   print(join("\t",
         $1,
         $2,
         $4?$4:($2),
         ($8 eq "ins")?$9:($5 eq "dup")?"-":$6,
         (($5 eq "dup") || ($8 eq "ins"))?$6:($9)?$9:"-")."\n")
}';
chr7    140453145       140453145       A       T
chr7    140453136       140453136       A       T
chr7    140453135       140453136       TT      CA
chr20   31022287        31022287        -       A
chr19   13054614        13054614        G       -

Regular Expression
This uses enclosing groups to describe the important bits of the sequence. The groups are numbered from 1 upwards, with each number corresponding to a left bracket (in order), capturing the entirety of the sub-expression until the paired right bracket. I'll describe the groups, which are referenced in the subsequent part of the script as $<group_number>:

Any number of characters, up to (but not including) the first ':'
Any number of the characters '0' to '9' [i.e. numbers]
Optionally, a "_", followed by any number of the characters '0' to '9'
The numbers from Group [3], i.e. excluding the "_"
Optionally, a "dup" or a "del"
Any number of characters selected from "A", "C", "G", "T" (i.e. bases)
Optionally, an "ins" or ">", followed by any number of bases
The "ins" or ">" from Group [7]
The bases from Group [7]

Print Code
The print statement prints out a tab-separated list, finishing off with a line break. I'll describe each list element:

Group [1]
Group [2]
If Group [4] exists, then Group [4], otherwise Group [2]
If Group [8] is "ins", then Group [9]; [otherwise]
If Group [5] is "dup", then "-", otherwise Group [6]
If Group [5] is "dup", or group [8] is "ins", then Group [6][otherwise]
If Group [9] exists, then Group [9], otherwise "-"

+1, although I'm tempted to quote Jurassic Park here: "You were so preoccupied with whether or not you could, you didn't stop to think if you should."

Here is an R solution, would probably be slower than the suggested Python solution though:
bed <- readLines("bed.txt")

res <- data.frame()

for(i in bed) {

  temp <- strsplit(i, ":") %>% unlist()

  # case N>N

  if(grepl(">", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer() - 1
    end <- begin + 1
    ref <- str_extract(temp[[2]], pattern = "[A-Z]")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del_ins

  if(grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- str_extract(temp[[2]], pattern = "_[:digit:]+")
    end <- substring(end, 2, nchar(end)) %>% as.integer()
    ref <- str_extract(temp[[2]], pattern = "[A-Z]+$")
    alt <- str_extract(temp[[2]], pattern = "[A-Z]+")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case dup

  if(grepl("dup", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- "-"
    alt <- str_extract(temp[[2]], pattern = "[A-Z]$")

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

  # case del only

  if(grepl("del", temp[[2]]) & !grepl("ins", temp[[2]])) {

    begin <- str_extract(temp[[2]], pattern = "^[:digit:]+") %>% as.integer()
    end <- begin
    ref <- str_extract(temp[[2]], pattern = "[A-Z]$")
    alt <- "-"

    temp_res <- data.frame(temp[[1]], begin, end, ref, alt)
    res <- rbind(res, temp_res)
  }

}

print(res)

  temp..1..     begin       end ref alt
1      chr7 140453144 140453145   A   T
2      chr7 140453135 140453136   A   T
3      chr7 140453135 140453136  TT  CA
4     chr20  31022287  31022287   -   A
5     chr19  13054614  13054614   G   -

Thank you all very much :).


text:
Find rsIDs for GRCh37 SNPs and rsIDs for GRCh38 SNPs and compare for overlap
I have a long list of variant IDs that I generated as a result of research done one GRCh37 genomes (e.g. 13_28025615_G_C_b37). I want to get their rsIDs and compare them to a newly-obtained list of GRCh38 rsIDs that have also been retrieved from similarly formatted variant IDs (e.g. chr22_50473581_G_A_b38). How would I go about doing this? I've tried installing biomaRt as a package but I can't quite make heads or tails of it and it's not clear to me that it would support both versions of genome builds.

it's not clear to me that it would support both versions of genome
  builds.

You can specify what version of the database you are querying.
I'd start with biomart on the ensembl website, see if you can get it to do what you want, then work on getting biomaRt to do what you want.

https://www.ensembl.org/biomart/martview/ > CHOOSE DATASET > Ensembl Variation > I can't find any option for GRCh37, only 38

The website might not have that enabled, but you can choose old versions in R.  snpmart = useEnsembl(biomart = "snp", dataset="hsapiens_snp", version = 90)

When I do that, the only datasets I can see with listDatasets(snpmart) are GRCh38.p10.

Okay, I had to go all the way back to version 75 from Feb 2014 to get GRCh37. Thanks for your help.

Don't use the e75 version of BioMart, it's really out of date. There's a dedicated GRCh37 version where the variation database has been updated. More info here: https://www.ensembl.org/info/data/biomart/biomart_r_package.html

rsIDs change with the dbSNP builds, not with the reference genome builds. rsIDs do not depend on a reference genome. They point to a specified locus regardless of the differences in genomic assemblies. rsIDs numbersare assigned by dbSNP database maintainers. The numbers may be deleted from database (e.g. when a submitter withdraw data) and may change on merges (e.g. when multiple rsIDs were assigned at the same genomic location).
Specifying the location by rsID is stable, and can be only affected by an eventual merge or deletion, if any, with a new dbSNP database build. It cannot be affected by a reference genome build.
However, a base position (offset) of a locus within a chromosome may change between the builds of the reference genome. So, when specifying the location by the chromosome and base position, the reference genome build should be also given, to locate an appropriate rsID.
You can download dbSNP database files (in VCF format) of the database build that was actual at a time when a specific build of a reference genome was also actual, but these files are huge. A possible place to download is https://ftp.ncbi.nih.gov/snp/pre_build152/organisms/ or https://ftp.ncbi.nih.gov/snp/archive/
However, these files are arranged by the dbSNP database build (e.g. 152, 153, 154), not by the reference genome build (e.g. 37 or 38).
If your list is relatively small, you can fetch the data rs numbers via dbSNP API by chromosome and location. The API is available at https://www.ncbi.nlm.nih.gov/snp/docs/eutils_help/
For instance, do “advanced search” of dbSNP at https://www.ncbi.nlm.nih.gov/snp/advanced
When it comes to search by the base position, you can only specify positions for GRCh37 (called “Base Position Previous”) and GRCh38 (“Base Position”). This URL will help you construct the request string that you will further use for the dbSNP API requests.
You may also need to add “NOT MergedRS[All Fields]” to exclude rsIDs that have already been merged. I wish I knew what to type instead of “All Fields” here to refine the search, but the documentation is very sparse.
In your examples, for 13_28025615_G_C_b37 specify Chromosome=13 and Base Position Previous=28025615 in the “Advanced Search”.
As about the API, for 13_28025615_G_C_b37 query https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=13[Chromosome]+AND+28025615[Base+Position+Previous]+NOT+MergedRS[All+Fields]
and for chr22_50473581_G_A_b38 query https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term=22[Chromosome]+AND+50473581[Base+Position]+NOT+MergedRS[All+Fields]
The first rsID that you will get will be the main current rsID, the others (if any) will be the merged ones.
Here is a small perl example on how to retrieve data from the dbSNP. The fetch_GRCh37 function returns the rsID for a given chromosome and base position. For GRCh38, remove the “+Previous” from the URL. You can call the fetch_GRCh37 in a script multiple times. It’s better that you get rsIDs for your both GRCh37 and GRCh38 data to make sure that you refer to the same rsIDs that are actual after eventual merges.
Please make sure to read and comply to the requirements before using the API.
A sample Perl script snippet follows:
use XML::LibXML;
use LWP::UserAgent;

sub fetch_GRCh37 {
        my $chromosome = $_[0];
        my $basepos = $_[1];
        my $url = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=snp&term='.$chromosome.'[Chromosome]+AND+'.$basepos.'[Base+Position+Previous]+NOT+MergedRS[All+Fields]';
        my $UserAgent = new LWP::UserAgent ();
        my $GET = $UserAgent->get($url)->content;
        my $dom = XML::LibXML->load_xml(string => $GET);
        my $numberofnodes = int($dom->findnodes('eSearchResult/Count')->to_literal());
        if ($numberofnodes > 0) {foreach my $title ($dom->findnodes('eSearchResult/IdList/Id')){return("rs".$title->to_literal());}}
}

print fetch_GRCh37(13, 28025615);


text:
How to map selected genes to Metabolic pathway Maps
I have a selected Arabidopsis Genome Initiative (AGI) list for RNA seq and proteomics data, how can I map them to metabolic pathway maps to vilualize(e.g. TCA cycle / FA / Photosynthesis) in KEGG or AraCyc or any freely available pathway maps?
Could someone share an example if possible?

KEGG has a webtool for this, KEGG Mapper.
You can just copy-paste your gene ids in there and see how your pathways are colored.

Thank you Benn :)


text:
What does final generated models with TM-score of DMPfold signify in "sequence.ema.txt" file?
So, my runs have finished and got 5 final models. I predicted TM-score of one model to understand the result. The output of "predict-tmscore.sh" for "final_1.pdb" is:
1 models read from PDB file
1 models read from PDB file
nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is -0.082424
Probability of each TM-score bin (cols 1-10), estimated TM-score (col 11) and prob TM > 0.5 (col 12) written to sequence.ema.txt
Estimated TM-score of ../final_1.pdb is 0.5

could you please help me understand what does "nres", "neff", "dpscore", and "dpscorenorm" below mean:
    nres is 377
neff is 15423
dpscore is -300.599030
dpscorenorm is    -0.082424

and it also generates a "sequence.ema.txt" file with below value:
2.0308916646172293e-05 2.0682286958617624e-06 0.04564645141363144 0.14857041835784912 0.27177685499191284 0.31814685463905334 0.19150513410568237 0.022412151098251343 0.001872337656095624 4.743059980683029e-05 0.5036169363265117 0.5339840054512024

Please let me know what these above values signify?
And,I would really appreciate your comments on what does "rawdistpred.1" file mean? There are total 9 such files. When predicitng TM-score, do I need to use "rawdistpred.1" file for all final models or I can use any "rawdistpred.*" file?
thanks!

Excellent question and answers on two core structure prediction programs

This is described briefly in the paper. The four quantities are inputs to the estimator of model accuracy:

nres is the sequence length.
neff is an effective sequence count, or how deep the alignment was.
dpscore is how well the final model matches the initial distance prediction.
dpscorenorm is a normalised version of the above.

Don't worry too much about directly interpreting these. They are only given for completion. The important line is Estimated TM-score of ../final_1.pdb is 0.5, indicating that DMPfold just about thinks this is a decent model.
The "sequence.ema.txt" is more raw output giving the probability of the TM-score being in each bin. "rawdistpred.1" is the distance predictions from the first iteration. Always use "rawdistpred.1" as input to the TM-score predictor as this is how it was trained.


text:
Protein model quality assessment
How do I perform protein model quality assessment?
I will obtain predicted structures from DMPfold, I-Tasser, and Rosetta, what should be the next approach to select a structure from among this set of models?
Note, I have been looking at some, appropriate methods such as MetaMQAP, QA-RecobineIt, but their servers are not working.

Please don't ask "what is the best". The best will always be a matter of opinion and the best method for your data will depend on exactly what kind of data you have. Please [edit] your question and rephrase it so it isn't asking for opinion (best is an opinion).

Edited to be compliant

Have a read of the Estimation of model accuracy in CASP13 paper to get the latest rankings.
Generally you would run the model quality software on the models you have generated and choose the best. Also have a look at the structures in PyMol, you learn a lot that way. See, for instance, how much variety there is, and where the models vary. Is there a conserved core in many of the models? Or is there just a blob of "spaghetti" that doesn't look like an ordered protein?
If you have many models (e.g. from a Rosetta run) you might want to look into clustering the models first to reduce the number you have to work with.

Co-variance with convolution is an extremely hot method, so DMPFold via psicov. Checking the reality of the folding is obviously important, convolutions can be missed, inaccuracies can happen.

DMPfold uses PSICOV in its pipeline. Does this mean that DMPfold does quality assessment in its pipeline?

It would be good for you to read the publications on psicov, to understand it. I believe @jgreener is part of the group that developed it

DMPfold has a model quality predictor that is specific to DMPfold. It can only score DMPfold models because it uses some internal measures, such as how well the final model matches the initial distance prediction. DMPfold shouldn't therefore be thought of as a model quality assessment method. PSICOV is part of DMPfold but isn't related to model quality assessment at all - it predicts residue-residue contacts by inverting covariance matrices.

Straight covariance results in convolution in the non- deep learning sense, ie pure statistical sense (I haven't explained that very well,but hope it makes some sense)


text:
Convert FPKM to RPKM
I have my 24 samples analyzed by Hisat2>Stringtie, from paired end sequenced data. I have got FPKM data and am planing to convert them to RPKM. to be able to compare my samples with another already published.
1- How can I convert from FPKM to RPKM?
2- Is it an acceptable method to compare samples using RPKM between two different studies or I need more or different procedures?

I don't know if you have followed the discussions from the last years, but FPKM and RPKM are obsolete. Statisticians found that FPKM and RPKM are not the right values that can be used for proper statistics. It is a wrong way of normalizing RNA-seq data. In stead they suggest to use raw read counts, which can be used in tools like edgeR or DESeq2. Also, CPM would be better than FPKM or RPKM, for example you can use (log2) CPM values in limma with trend.
But to answer your questions.

You can't convert, you need to recalculate from raw data.
No not acceptable. To compare studies you need to correct for batch, which you can do with statistical programs such as edgeR and DESeq2.

The best advice I can give is to start with raw data, also from the already published data that you would like to compare your own experiment with. And continue with edgeR or DESeq2 or the like.


text:
Metagenomic shotgun data with internal control
For 24 human stool samples I have metagenomic shotgun reads from an Illumina platform. The reads were filtered and mapped against a bacterial species library and specific maps to species were kept and counted. Importantly, before DNA extraction from any sample it was spiked with a solution of a fixed DNA concentration of Thermus thermophilus, a species not occurring in the human microbiota. The reads from T. thermophilus allow me to convert relative abundance estimates back to abundance in the stool sample, to absolute abundance.
I want to test whether T. thermophilus is indeed a good internal control that removes confounding by converting to absolute abundance. To test this, I compare between-species correlations on the relative abundance scale with their correlations on absolute abundance scale. I would expect less correlation in the latter. I first take the raw read counts of a random set of species, I divide counts by a sample-specific size factor estimate (total number of reads r in sample i / geometric mean of r's) to adjust for library size, then I do arcsinh transform (to stabilise variance) and then I look at between species correlation. Next I take again the raw read counts, correct by the size factor, then correct for an estimate of the coverage ( estimated proportion of T. thermophilus in sample i ), and I do the acsinh transform. Again I look at between-species correlations.
I compare species-pair correlation before and after T. thermophilus correction in a scatterplot (see attached) to visualise the effect of the correction. I expected that the point cloud indicates a slope <1 reflecting that previously strong correlations have weakened on the absolute scale. But that is not what my result shows. Indeed, some weakly positive correlations have become stronger. How does this make sense? What other ways are there to test validity of an internal control?

The problem with correlations in relative scale is know on the literature, it has been suggested to change to other metrics. So I wouldn't use it for internal control.
As you can see on figure 3 of the linked article the correlation on the absolute scale tend to be higher. So I would say this is within expectations.
As other way for internal controls are blanks (Just water to see if there are any contaminants), and technical replications, to make sure that the measure is correct not some problem with picking some amplifications instead of others...
Also you don't mention it but the tool and library used to impute the microorganism is also important. You can check with other tools to see if the estimation is similar.


text:
Best way to "fish" long reads that match a query sequence
Very simple set-up: We have a PacBio long high-fidelity (HiFi) reads, genome sequencing, and some of those that contain a particular sequence. I want to find out which ones. The length of the query sequence is 7kb, though ideally I would use as query a sub-sequence of only 1.5kb.
The output can the just the read ID or the full read sequence.
We tried DAMAPPER and it gives some decent results but it might not be appropriate to search for the shorter sub-sequence of 1.5kb. It also feels like an overkill for this purpose.
Bear in mind that this is my first foray into long-read sequencing so my mindset is still very much one of short, low error-rate reads.

I have seen the blasr aligner typically used for aligning PacBio reads to a reference.  Is there a reason why blasr would not be a good choice? I guess the very short sequence you are trying to align to might be an issue.
Minimap2 is capable of mapping pacbio reads, and can also do pairwise mapping to see if reads overlap with each other, so this function might help you:
From https://github.com/lh3/minimap2
minimap2 -x ava-pb  reads.fq query_seq.fq > ovlp.paf


text:
bedtools intersect on very large .bam file - -sorted confusion
I need to do a bedtools intersect operation on a very large .bam file.  When I use the standard bedtools intersect operation, the process consumes all the memory on my system.
From the bedtools documentation page, I find the -sorted option to bedtools intersect, which sounds like it is the perfect solution to the problem, but the documentation seems incomplete.  From the documentation, I need to ensure that the .bed file is sorted and I can use either bedtools sort or the Linux sort command to accomplish this:
bedtools sort -chrThenSizeA -i input.bed > sorted.bed

sort -k1,1 -k2,2n input.bed > sorted.bed

If I then execute the intersect operation on the .bam file, I get an error.
bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out
ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files.

Tracking this down I find that samtools sort, which was executed on the .bam file, "sorts by the leftmost coordinate, or the read name if -n is used."  It would appear there is no way to put the .bam file in order by reference name, which is what bedtools instructs.
For a workaround, bedtools sort does have the -g option, which allows me to specify the ordering of the .bed file.  To get this ordering, I execute:
samtools view input.bam | cut -f3 | unique > order.txt
bedtools sort -g order.txt > alt_sorted.bed

This puts the .bed file into the same order as the .bam file, and the bedtools intersect -sorted command now works, but is it working correctly?  It seems that bedtools intersect -sorted wants the files in consistent ordering, which can be done, but is that consistent with the intersect algorithm itself?  From the documentation's instructions to sort by the first 2 fields, chromosome name and left-most coordinate, it seems inconsistent to order it as the .bam file is ordered, by read name and not reference name.
I'm not sure how to check this assumption.
EDIT:  Some extra details regarding sorting of the .bam file.
If the .bam file is sorted with samtools sort input.bam > sorted.bam I get the error, 'ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files`
If the .bam file is sorted with samtools sort -n input.bam > sorted.bam I get this error, Error: Sorted input specified, but the file sorted.bam has the following out of order
It would seem that if I sort without -n, it thinks the chromosomes are out of order, and if I sort with -n option it thinks the coordinate order is wrong. The only fix I have is forcing the order of the .bed file with the -g option, but again, I have no idea if that is consistent with the underlying algorithm's assumptions for the -sorted option.

As a workaround, you could just split your bed file into one file per chromosome, run the intersect for each of those, and then merge the resulting bam files. Clearly not optimal, but if all else fails, it might give you a way out.

Thanks for the suggestion.  Unfortunately, that's not a workable option given the workflow that I currently have.

It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the -n option in samtools sort. Otherwise, if you used samtools sort without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered).
I think you are getting confused about 'ordering' the reference names. There are multiple ways to sort names; Some people might want numeric sorting (chr1, chr2, chr3...) and others might want lexicographic sorting (chr1, chr10, chr11...). Because of this, many tools allow you to specify the way you want to order the reference names.
There isn't necessary only one correct way to sort the reference names; what matters is that the order is the same between the files you're trying to intersect. This is exactly the purpose of bedtools intersect -g.
If you want to check the order of reference names in your bam, you can easily do so by checking the header via samtools view -h input.bam. It will say something like this:
@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1     LN:248956422
@SQ     SN:chr10    LN:133797422
@SQ     SN:chr11    LN:135086622
@SQ     SN:chr12    LN:133275309

The @SQ lines tell you how your bam is sorted. My suspicion is that this order is different from the Linux sort or bedtools sort default, which is why your bedtools intersect was failing.

It would seem that if I sort without -n, it thinks the chromosomes are out of order, and if I sort with -n option it thinks the coordinate order is wrong.  The only fix I have is forcing the order of the .bed file with the -g option.

samtools sort -n does not sort by reference name (the 3rd column of a bam file). It sorts by read name (ie. the first column of a bam file).

Inadvertently deleted my first comment. while updating my original post.  Oops.   Here's what was in my clipboard buffer: Yes, I've done sorting with samtools sort without the -n option, giving the error above. If I sort with the -n option, I get a different error: Error: Sorted input specified, but the file sorted.bam has the following out of order If I look at the header of the .bam file when the -n option wasn't used, the first line is '@HD VN:1.0 SO:coordinate` . I interpret SO:coordinate to mean that the Sort Order is by coordinate, which is consistent with what the documentation says.

As per Lynn's answer, sorting is dependent on how chromosomes/contigs are named, how you want them to be sorted, etc.
samtools sort or sambamba sort will sort as per the SAM specification by the order of the contigs in the header of the file. If you've aligned these sequences to a human reference genome, typically it will be chr1, chr2, ... chr10, chr11, ...
To sort the BED file in the same way, you can specify the sort order with the bedtools sort -g option and providing the explicit order.
Alternatively, use sort with the -V, --version-sort option, that provides a "natural sort of (version) numbers within text".
sambamba sort -o output.bam input.bam
sort -k1,1 -V -k2,2n input.bed > sorted.bed
bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out

Thanks for the response.  As I mention above, I have used the -g option that you suggest, which seems the most useful solution for my workflow.  The concern I had is that this specified ordering is not consistent with how bedtools would normally sort the .bed file and question is whether this atypical sorting would affect the intersect operation.

As far as I'm aware, bedtools doesn't care what positional order you sort them in, just that all the files are in the same order. So bedtools sort -g is probably your best option to force a specific chromosomal ordering

Thanks, James!  That is consistent with what I see, but not knowing what's going on under the hood makes me question things a bit.


text:
Getting species from ensembl sequence id
I have downloaded some orthologues from ensembl website in fasta format. The problem is that it uses the ensembl sequence id as name for the sequences. I want to know the species name or any other information about the sequence (chromosome, location...) For example for this sequence id: ENSMODP00000018933. How can I do that? I have been looking at the REST API but the most I can find is an API call for getting the sequence passing the id as parameter.
https://rest.ensembl.org/documentation/info/sequence_id
I have also tried this API call but does not work for sequence ids:
https://rest.ensembl.org/documentation/info/xref_id

Use efetch for NCBI

The lookup/id endpoint will get it for you. Of you can just look up the three letter species code, in this case MOD, on the list in the documentation.

I made a script that uses the list in the documentation. The lookup returns me almost always that: {"error":"Expand option only available for Genes and Transcripts"} Thanks.


text:
awk working with large files
I have two very large vcf files 2GB and 6GB I want to look for unique combinations of CHROM and POS and output the row that matches. However, because the files a so large my machine always hang and stop processing. Is there are a way to work around the problem of this large file?
I am using this command, taken from an answer to another question:
awk  '{ if(NR==FNR){a[$1$2]=$0}else{if($1$2 in a){print}}}' file1.vcf file2.vcf

file 1
#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   10611   rs189107123 C   G
1   51954   rs185832753 G   C
1   13327   rs144762171 G   C

file 2
#CHROM  POS ID  REF ALT QUAL
1   10366   rs58108140  G   A
1   51935   rs181754315 C   T
1   51954   rs185832753 G   C
1   52058   rs62637813  G   C
1   52144   rs190291950 T   A

output
1 10366 rs58108140  G   A
1 51954  rs185832753 G   C

Do you have actual vcf files? The lines you show are not valid vcf, is that really what you have or do you have proper vcf files?

Maybe comm, which is part of GNU coreutils, is more efficient.

comm writes to standard output lines that are common, and lines that are unique, to two input files; a file name of ‘-’ means standard input. Synopsis:
comm [option]… file1 file2

In your case you could sort both files based on the first 3 columns:
file1=file1.vcf
file2=file2.vcf

comm -12 <(awk '{print $1,$2,$3}' $file1 | sort) \
         <(awk '{print $1,$2,$3}' $file2 | sort) > matches

# result:
cat matches
1 10366 rs58108140
1 51954 rs185832753
#CHROM POS ID

This could be of course polished by removing the header.
Also if you want the original lines as output, you can grep the matches to the original files:
grep -f matches file1
#CHROM POS ID REF ALT QUAL
1 10366 rs58108140 G A
1 51954 rs185832753 G C

Edit: In the first version, intermediate temp files were used, but as suggested by terdon in the comments, this can be avoided using Process Substitution allowing the <(..) contstruct.

You can also avoid creating temp files with comm -12 <(awk '{print $1,$2,$3}' file1 | sort) <(awk '{print $1,$2,$3}' file2 | sort).

That's true, Process Substitution avoids temp files very nicely, but since OP was talking about large files this might not fit in the memory, so writing to disk might be more appropriate for that particular scenario.

I doubt that would be an issue. sort will create temp files when sorting anyway, this won't be a memory intensive procedure.

If that's the case, I'd add it to the answer. Do you know what the official procedure is for that? Is it ok to include this as on top as the quickest solution, or should we just leave it as your comment here?

Oh, it's absolutely OK to include anything given in the comments! More than OK, it's actively encouraged. The only reason I didn't edit it in myself was that I wanted to let you choose whether or not to include it.

Good to know! Hope this is ok this way now! Thanks

Sort both files by chromosome and position and read each file one line by line.
Either you iterate to the next line in one file, or go the next in the other, or there's a match that you print out and then iterate both lines.
This is not necessarily an awk approach, but is doable with Perl, Python, or any other scripting language.
In pseudocode:
Sub getCommonLines()
  Sort FileA
  Sort FileB
  Get lineA from FileA
  Get lineB from FileB
  While not EOF in FileA or EOF in FileB:
    If lineA > lineB:
      Get lineB from FileB
    Elseif lineB > lineA:
      Get lineA from FileA
    Else:
      Print lineA or lineB
      Get lineA from FileA
      Get lineB from FileB
End Sub

The use of the > operator is based on lexicographical and numerical ordering provided from sorting files A and B, and comparison should be done only on the fields of interest. If you compare other fields (esp. from VCF files) you may get other fields that give you the wrong or unexpected comparison result.
The term EOF means end-of-file, which means just what it says — you're at the end of the file and there are no more lines to compare between files, so your subroutine is done.
Sorting has a time cost, but memory overhead for modern sorting algorithms is very minimal. You can sort via process substitution, but sorting to an intermediate file can be useful if you have other things you want to do with them, where sorted input is valuable.
This would be a good approach for very large files that you do not edit or update often (or at all), for which a hash table (or "dictionary", or "associative array", all the same terms for the same thing) is impractical because of its large memory overhead.
You can also use a tool like comm with sorted files, as described in another answer/comment. This does basically what my answer describes without you having to understand or implement things. You must remember to provide sorted input, however, as with other Unix tools that do line comparison (uniq, etc.).
Despite this answer being a really long "comment", it is useful to have a basic understanding of how things work under the hood, so that when you use these tools, you know what's going on and why you need to provide inputs in a certain way. It's easy to get a million one-liners off of Stack Exchange sites, but understanding is valuable beyond gold.

Can you expand this to show how it could be done? As it stands, this is really more of a comment than an answer.

Use bedtools intersect. You can select the various combinations of outputs:
https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html


text:
Are there databases or resources that contain expression of proteins or genes before and after cancer treatment?
Are there databases that contain information about gene or protein expression before and after targeted therapy or chemotherapy?
I am curious about different cancer cell lines or cancer types.
Any suggestion on a set of resources I should look into?

"Are there databases..." As a general rule, no, there are almost never specialized databases for anything.

@DevonRyan, any idea how to extract expression ratio before and after cancer treatment from https://portal.gdc.cancer.gov/, https://bioinformatics.mdanderson.org/public-datasets/, https://www.immport.org/shared/home, https://cancer.sanger.ac.uk/cosmic, or http://www.cbioportal.org/

Find appropriate samples, match them up, find expression metrics, write a program to do that. Never expect tools for random things like this.

One such database/tool is “[PharmacoGx][1]: an R package for analysis of large pharmacogenomic datasets.” Bioinformatics (Oxford, England).
They brought together a large compendium of experimental data (cell lines and the like if I remember correctly from a talk), standardized these and made available within a special object class. I think there is genomic and transcriptomic data, not sure if there is proteomic data though.

Do you mind to include a simple example?

What you are trying is not a simple task overall, therefore I cannot provide a "simple" example. I can only point you to their vignette, they have worked-out cases there: http://bioconductor.org/packages/release/bioc/vignettes/PharmacoGx/inst/doc/PharmacoGx.pdf

I see. What about the databases COSMICS, TCGA, and immport. Aren’t they supposed to have such things?

Can't talk about others but TCGA is (mostly) treatment naive, you will need to check clinical information tables (of your cancer of interest) if patients got treatment or not though. If you decide go down that road, TCGABiolinks (and especially its vignettes) would be very helpful.


text:
How to edit the headers of multiple fasta files from multiple folders?
My directories are organized as follows: one main directory, in which I have multiple directories that end with a number ranging from 314 to 727, followed by .3   . For example, 'mgm4761314.3'. Within these directories, I have fasta files named after a similar pattern (for example, mgm4761314.3.299.screen.passed.fna is in the directory mgm4761314.3).
There is only one such file (.fna) per directory.
What I've been trying to do is to write a loop that would go through all these folders and fasta files, and add the short number associated with the folder/file name (eg. 314, 315...727) to the headers. For example, for my first folder/fasta file, instead of:
>seq1
>seq2

I want:
>314_seq1
>314_seq2

The output would be exported to a file called, for instance, 314.fna.
I've sort of made it work for one number at a time, but I have plenty of files, so I'd appreciate if anyone could help me with the loop part! This is what I tried (and it didn't work). Using ubuntu.
for i in {314..727}; do sed "s/>/>${i}_/" [...]/*${i}.3/mgm4761${i}.3.299.screen.passed.fna > [...]/*${i}.3/${i}.fna; done

This should be easy enough, but please [edit] your question to clarify the directory structure. Does this need to be done to all files with an .fna extension in each directory?

@terdon Edited! Basically, there is only one file per directory (and all files are .fna)

Or could go perl -p -i -e 's/

You were so close! This should do the trick:
for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done

Note that I am using -i which makes sed edit the original file. Because I use -i.bak, this will create a backup file with the original data. Once you are sure the command does what you need, you can delete all of these files with rm */*fna.bak.

This seems to be a step in the right direction! However, there's something about it which went weird. It does what it's supposed to, but it also duplicates the headers. Example: '>314_NS500799:83:H53J2AFXX:2:11101:11350:1048

314_NS500799:83:H53J2AFXX:2:11101:11350:1048
GGGCAGTCCTTGCAGGGTGATACAGTCACCGTTTCTCTGAATGGGGAGGATCGTCCCCAAACTGTCAAGTACTACCGTTTGGCCCAGATCGTCTTTCAGAACCCGTACTCTTCACTGAATCCGCGCAAGCCGGTGCGCGACATCGTCGGCA
314_NS500799:83:H53J2AFXX:2:11101:2772:1050
314_NS500799:83:H53J2AFXX:2:11101:2772:1050
etc'

IMO perl -p -i -e gives you more regex power, but this problem is easily solved just anchor the match and delete to the end

@Laura argh! That's because I'm an idiot and left a debugging command I had used while testing my answer. If you've run it already, you will want to rename all the bak files back to the original. You can do this with: for i in */*bak; do mv "$i" "${i//.bak}"; done. Then, run the updated command (without the p at the end of the s///) and it should work as expected.

@MichaelG. yes, perl is indeed better (and is where sed got the -i from), but overkill here. The match is already anchored to the beginning. The issue was that I had used s/foo/bar/p and the p causes sed to print any line where the substitution was successful. I was using that along with -n which suppresses output to test that it was working as expected, but forgot to remove it when posting the answer.

@Laura if you didn't keep the backup files, you can remove the duplicated headers with perl -i -ne 'if(/^>/ && $_ eq $last){next}; print; $last=$_; ' */*fna.

@MichaelG. Do you know, by any chance, how I could add a layer of complexity to this code? As in, instead of adding from a range of numbers that happen to be equivalent to folder endings, would there be a way to do this if I didn't know those numbers from before? I tried this to add codes from multiple folder names to the fasta, but...it's so bad that I'm not even getting an error message.
for j in ${j}; do         for i in ${i}; do                 sed -i.bak "s/^>/>${j}_${i}_/" mgp${j}/mgm????${i}.3/*.fna         done done

@terdon is the expert on bash and awk. Personally I'd write a Perl script using glob, you might achieve it by piping.

Looking at your code, that is  alot  of directories why not replace ??? With *

@Laura I'm afraid I can't answer without knowing your directory names. But as Micahel said, for dir in * might be enough. You might want to ask a new question though where you can give us the relevant details (the devil's in the details with this sort of thing).


text:
How to map multiple CDS sequences to exome data using bwa mem?
Normally, I use bwa-mem to map my fastq reads for a bunch of species using one protein-coding gene CDS from my reference species. This was fine when I had to try it for just a group of genes...I would obtain single CDS fastq files and use them as a reference.
However, we now want to obtain protein-coding sequences for ALL GENES in my exome datasets. If I have one file containing all the CDS from my reference species, can that be used in bwa-mem? The end of my custom script produces one multiple sequence alignment file when I ran it using one gene only...my concern is that bwa-mem will somehow replace certain CDS or would I just end up with over 20000+ MSAs?  Which is fine because I need them to be individual MSAs for each protein-coding gene.
I use a simple python script at the end to extract base information from the mpileup files generated to create my MSAs.

If I have one file containing all the CDS from my reference species, can that be used in bwa-mem?

That's fine, normal usage of bwa-mem entails giving it a fasta file with all chromosomes/contigs/transcript/etc. in an organism, so your usage is no different. How your custom "make an MSA" script handles that is unknown of course.

Yes unfortunately I need to figure out a new python script - I had tested it out with a small reference file containing about 10 CDS, but my final MSA was messed up it looked like it tried to merge all the CDS into one sequence. It made no sense at all. And my reference file contained a list of CDS so it was a list of 10 individual fasta sequences of diff lengths, not one fasta sequence concantenated. Maybe that was the problem? I am not sure how to make bwa-mem read each CDS line as its own reference sequence.

I am trying a new way out now where I assemble my reads against each chromosome from my ref species...but then I have no idea how to pull out all the individual protein-coding gene sequences. I am having trouble thinking of this differently since I have been running individual genes for the past year with no issue.


text:
How to obtain the extended unaligned bases of the read in reference based genome assembly?
I did illumina read mapping based on a particular gene sequence (reference sequence) using the following command,
bowtie2-build -f ref_gene.fasta ref_index && bowtie2 -p 4 -N 1 -t -x ref_index -1 genome_R1.fastq -2 genome_R2.fastq -S tal.sam && samtools view -bS tal.sam | samtools sort - -o tal.bam && samtools mpileup -uf ref_gene.fasta tal.bam | bcftools call -c | vcfutils.pl vcf2fq > cns.fastq && seqtk seq -aQ64 -q20 -n N cns.fastq > contig.fasta

It maps the reads according to the reference sequence and yields, exactly same length of aligned/assembled sequences. But, I need obtain the unaligned overlapping read bases at the both end (reverse and forward end)  Please see the image for better understanding.

Is there any way to get the same?

You're not going to get vastly different replies here than Biostars or github.

@DevonRyan sniper!


text:
How to convert a list of genes to one element per row?
I got two lists of about 300 genes from my GO analysis, like the following:
Mtss1/Sox11/Fbn2/Mmp2/Tcf7l2/Six2/Ahr/Fat4/Lef1/Stc1/Acvr2a/Glis1/Nog/Igsf10/Twist2/Igf1/...

How can I convert the list to one element per row (for Venn diagram  on Venny 2.1), as:
Mtss1
Sox11
Fbn2
Mmp2
...

$ echo "Mtss1/Sox11/Fbn2/Mmp2/..." | tr '/' '\n' > list.txt

For more information on tr: https://en.wikipedia.org/wiki/Tr_(Unix)

It works well in my terminal, thank you so much!

In R you can use tidyr's function separate_rows.
install.packages("tidyr") ## if not already installed
library(tidyr)
GO_results<-read.table("genes.txt")
separate_rows(GO_results, GO_results$V1, sep = "/")

Many ways, assuming you're on a *nix system.

tr. Alex already gave you a fine tr solution, but here's another approach assuming you have an input file:
tr / '\n' < golist.txt > newlist.txt

sed
sed 's|/|\n|g' golist.txt > newlist.txt

Or, to edit the original file:
sed -i 's|/|\n|g' golist.txt

Perl
perl -pe 's|/|\n|g' golist.txt > newlist.txt

Or, to edit the original file:
perl -i -pe 's|/|\n|g' golist.txt

grep (GNU grep, the default on Linux systems)
grep -Eo '[^/]+' golist.txt > newlist.txt

Thank you all very much!

@Lin you're very welcome. If one of the answers solved your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. Since Alex answered first, you might want to accept his.


text:
Weird phylogenetic tree
I am required to create a phylogenetic tree(s) of my genes of interest and am currently blasting (tblastn) my genes of interest from S. mediterranea against model systems before expanding my search.
Problem I get something like the attached which does not look like right. What am I doing wrong? I'm aligning via MUSCLE, and running Maximum Likliehood Protein models and consequentially, using the correct model (JTT, for instance).
Greatly appreciate any help in advance.

Which gene is it?

Anyway yes it's possible, I'm confused by your Tblastn because that will give nucleotides, JTT is s protein model

@MichaelG. tblastn gives protein, not nucleotide. It takes a protein query and blasts it against a translated nucleotide database. The results are protein alignments.

@zaci, could you explain what looks wrong? Why would there be something wrong here? You seem to be aligning small parts of each protein, so it makes sense that you're not finding the "real" species tree.

I suspect it's gene duplication. Helminths and insects are in the right place as basal metazoa

Anyway the chances are this is gene duplication. It will be a housekeeping type gene or else histone that sort of thing. You simply map the different genes onto the various metazoan classes. In your case nematodes and arthropods are basal metazoan, so their is correct given the clade they form part of all all representatives is that gene.
There is nothing weird you just need to read the deep phylogeneticists work like Tom Caveliar-Smith and Sandy Balfour, possibly Blair-Hedges
There are alot of eukaryotes prior to the invertebrates, just look at the Cambrian explosion

How about the tri-furcation at the branch giving rise to RFX5-RFX7, Smed  RFX15486 and Dmel Uncharacterized...? Isn't there a problem there?

No, its simply a polytomy, increase the amount of data and resolution will sharpen to a bifurcation. From recollection the arthropod - worm boundary isn't well resolved. Might have changed with genome availability

I think Blair-Hedges was the person most related to this


text:
awk find Sickle cell anemia from father, mother, son, daughter
I have a file where i what to check for a pattern below and print out all the lines that matches the pattern but the output is empty.
here is the pattern that i am looking for.

Father = 1/0 or 0/1 or 1|0 or 0|1
Mother = 1/0 or 0/1 or 1|0 or 0|1
Daughter 1 = 1|0 or 0|1 or 1/0 or 0/1
Daughter 2 = 1|1 or 1/1
Daughter 3 = 0|0 or 0/0
Son 1 = 1|0 or 0|1 or 1/0 or 0/1
Son 2 = 1|1 or 1/1

mysample file
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9
1 762005 . G C 1/0:36:1.191:5.321:0.139:360.9:357:0.00:1.93:0.97:0.00:C,18,2.166,G,18,3.156:18,18 0/0:24:0.799:4.207:0.167:0.0:60:0.00:0.36:3.26:0.00:G,24,4.207:24,0 1|0:40:1.332:3.384:0.075:332.2:332:0.22:13.83:0.22:0.00:C,19,1.186,G,21,2.197:21,19 1|0:31:1.032:3.290:0.097:236.2:236:1.75:8.18:0.07:0.00:C,13,0.134,G,18,3.156:18,13 1|0:32:1.066:4.290:0.125:187.0:187:5.67:1.78:9.77:0.00:C,11,0.114,G,20,3.177,T,1,1.000:20,11 1|0:43:1.423:4.405:0.093:181.6:182:10.13:1.24:0.45:0.00:C,14,2.124,G,28,2.270,T,1,0.010:28,14 0|0:34:1.125:2.332:0.059:0.0:96:0.00:4.09:0.26:0.00:G,34,2.332:34,0

I tried this, but the output is empty:
awk  '$9=="1|1"  || $10=="0|0" || $12=="1|1" || $8=="1|0" || $11=="1|0" || $6=="1/0"||$7 =="1/0" {print}' variants.vcf > test.txt

Should it match all conditions or any conditions? If Daughter 2 = 1|1 or 1/1, should that line be printed irrespective of the other genotypes?

Yes should match all

So, print all lines where at least one of each of your bullet points is true? Then why are you using || in your awk?

The fields aren't only the genotype. They contain much more information. So $9 will never be exactly 1|1, what you are interested in is that the field starts with 1|1. Also, if you need all of the conditions to match (meaning that at least one of the criteria for each sample is met), you need && not ||. So try this instead:
awk '($6~/^1\/0/   || $6~/^0\/1/ || $6~/^1\\|0/   || $6~/^0\\|1/) &&
     ($7~/^1\/0/   || $7~/^0\/1/ || $7~/^1\\|0/   || $7~/^0\\|1/  ) &&
     ($8~/^1\/0/   || $8~/^0\/1/ || $8~/^1\\|0/   || $8~/^0\\|1/) &&
     ($9~/^1\\|1/  || $9~/^1\/1/ ) &&
     ($10~/^0\\|0/ || $10~/^0\/0/ ) &&
     ($11~/^1\/0/  || $11~/^0\/1/ || $11~/^1\\|0/ || $11~/^0\\|1/) &&
     ($12~/^1\\|1/ || $12 ~/^1\/1/)' foo.vcf

Running that on your file returns:
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9
1 762005 . G C 1/0:36:1.191:5.321:0.139:360.9:357:0.00:1.93:0.97:0.00:C,18,2.166,G,18,3.156:18,18 0/0:24:0.799:4.207:0.167:0.0:60:0.00:0.36:3.26:0.00:G,24,4.207:24,0 1|0:40:1.332:3.384:0.075:332.2:332:0.22:13.83:0.22:0.00:C,19,1.186,G,21,2.197:21,19 1|0:31:1.032:3.290:0.097:236.2:236:1.75:8.18:0.07:0.00:C,13,0.134,G,18,3.156:18,13 1|0:32:1.066:4.290:0.125:187.0:187:5.67:1.78:9.77:0.00:C,11,0.114,G,20,3.177,T,1,1.000:20,11 1|0:43:1.423:4.405:0.093:181.6:182:10.13:1.24:0.45:0.00:C,14,2.124,G,28,2.270,T,1,0.010:28,14 0|0:34:1.125:2.332:0.059:0.0:96:0.00:4.09:0.26:0.00:G,34,2.332:34,0

Hi @terdon, is there a particular reason for escaping the | with two backslahes?

@haci yes, the double escape is needed inside // in awk. Something about how awk interprets regexes in its match (//) operator. An unescaped | means "OR" in regular expressions, so it needs to be escaped. The double escape is an awk thing.

Could you please take a look at this question? https://bioinformatics.stackexchange.com/questions/10548/awk-find-pattern The answer I proposed for that question (which basically is more or less your answer above except for single backslashes to escape |) gives the expected output, whereas the one with double backslashes matches lines/columns it should not and I cannot figure out why.

@haci um. I don't know. I can't figure it out either. I know awk sometimes requires the double escaping, but I never remember exactly when or why. I also know that in this case, it doesn't work without the double escaping, but you're absolutely right, in your answer it fails with the double. I think I'll probably have to post a question on [linux.se], I just need to make a minimal reproducible example of this. I'll let you know, if I do.


text:
Filter out other character after the genotype
Is it possible to filter out all the characters that come after the genotype, for example, I have a data
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A 0/0:37:1.224:2.367:0.054:0.0:108:0.00:4.75:0.53:0.00:T,37,2.367:37,0 1/0:33:1.099:3.320:0.091:358.5:359:0.59:4.34:0.59:0.08:A,15,1.149,T,18,2.171:18,15 0|0:30:0.999:1.314:0.033:0.0:87:0.00:1.16:0.00:0.00:T,30,1.314:30,0 0|1:31:1.032:2.299:0.065:211.5:212:2.61:0.48:0.07:0.43:A,12,0.124,C,1,0.010,T,18,2.165:18,12 0|1:37:1.232:2.374:0.054:263.2:263:1.47:8.69:0.06:0.05:A,16,1.160,T,21,1.214:21,16 0|0:43:1.423:2.447:0.047:0.0:98:0.20:0.83:0.45:0.00:A,1,0.010,T,42,2.437:42,1 0|1:40:1.323:3.393:0.075:280.5:281:1.95:2.36:3.47:0.27:A,17,1.175,T,23,2.219:23,17
1 758663 . CGCC TTAG 0/0:27:0.893:0.296:0.000:0.0:90:5.15:0.85:3.67:0.01:27,0 0/1:21:0.699:0.234:0.000:2333.1:353:0.93:2.17:2.36:0.00:9,12 0|0:35:1.166:0.370:0.000:0.0:108:3.97:3.43:4.75:0.00:34,0 0|1:33:1.099:0.375:0.000:2820.8:447:0.07:8.69:0.53:0.00:17,15 0|1:23:0.766:0.232:0.000:1860.0:382:0.85:0.17:0.09:0.00:13,10 0|0:23:0.761:0.255:0.000:0.0:75:1.51:2.59:0.78:0.00:22,0 0|1:18:0.596:0.215:0.000:1686.8:275:0.00:2.17:0.10:0.00:9,9

i want to filter out everything that comes after the genotype and have an output like below
#CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2
1 736689 . T A         0/0     1/0     0|0     0|1        0|1     0|0   0|1
1 758663 . CGCC TTAG   0/0     0/1     0|0     0|1        0|1     0|0   0|1

Using sed:
sed 's,\([01][/|][01]\):[^ ]*,\1,g' file
#    ^  ^               ^     ^  ^
#    |  |               |     |  *--- applies to all found pattern
#    |  |               |     *------ print the content inside the parenthesis
#    |  |               *------------ any character except a space
#    |  *---------------------------- the pattern to be kept
#    *------------------------------- substitution

You also may want to format the output using column command:
sed 's,\([01][/|][01]\):[^ ]*,\1,g' file | column -t
#CHROM  POS     ID  REF   ALT   FATHER  MOTHER  DAUGHTER1  DAUGHTER2  DAUGHTER3   SON1  SON2
1       736689  .   T     A     0/0     1/0     0|0        0|1        0|1        0|0   0|1
1       758663  .   CGCC  TTAG  0/0     0/1     0|0        0|1        0|1        0|0   0|1

Note the comma separators in the sed command can replaced with other characters typically "/"


text:
Sequencing center returns data in several files
We send some samples to sequence and we got several (fastq.gz) files for each sample. The files are distributed at two or three folders with different dates (more than a week apart). The dates of the folders where the data is in:
2019-08-12
2019-08-21
2019-08-21
2019-09-04
2019-09-14

I asked the facility about why there are the same file in different folders, and they told us that:

we have split sequencing of that project on different flowcells, so all
  files are correct and you have to merge the individual (duplicate) files
  to get the requested total read numbers.

How can different flow cells be from the same experiment and be run with more than a week of difference? (I thought that all the flow cells are run at the same time in each sequencing)
Can be there some batch effect? How can I check if there is batch effect if we merge the files before trimming, mapping and counting?

Usually, a single sequencing library is created from each sample. Assuming you are using barcoded indexes to identify multiple samples, the index is attached to each library, and the libraries are pooled to make a library pool.  A library pool can be frozen, and is very stable.
A pool can be sequenced multiple times by pipetting a portion of the pool and supplying it to the sequencing process.  Each time, a sample of the molecules in the pool will be sequenced.
Generally, it is true that sequencing runs are very repeatable, and that you will not suffer any batch effects.
If you want to check for batch effects you would need to analyse all of your sequence files separately, and check this first, or in parallel, to your main analysis.
One thing you might want to check is whether the same library pool was sequenced multiple times, or whether separate pools or libraries were made for each sequencing run.  In the second case, there may be more likely to be batch effects.

Oh, thanks. I think they submitted a portion of the library pool because the index are the same. However, why they didn't add more sample from the pool from the beginning?

Regarding your sample pooling question, probably this is not happening in your case but here is our workflow, which generates multiple output files per sample like you get: We use 10x and the most expensive part is the sequencing run. As single cells are delicate and not all "droplet formation" steps go smoothly, first we do "shallow sequencing" on NextSeq to see if there are enough number of good quality cells in a given sample (after mapping, read counting, single cell QC) and then continue with deeper sequencing (one or more runs) with NovaSeq until we reach saturation with number of reads.

The number of sequences has a limit per sequencing run, even if you add more of the pool to the run - it is limited by the flow cell.  For this reason, it is normal to only sequence a part of the pool in any one run - enough to fill the flow cell - and if you need more reads, then carry out more runs with additional flow cells.

Another alternative would be to allocate some samples to one run, some samples to another run, and to fill up the flowcells that way.  However, this would mean that, if there are any batch effects it would be much more difficult to separate them from sample effects.  At least if all samples are included on every flow cell, the batch effects will not affect any between-samples comparison.

@Jonathan Well, in our case not all samples where included in all the runs, that's partially why I ask about the batch effect.


text:
How do I Investigate and test the hypothesis on the effect of street-light regime on insects?
I need advice on a suitable analysis of a dataset that has an akward design because we investigate a real-world situation. We  have two sites (A, B) where we can conduct the research. The hypothesis is: under-dimmed light, we catch less number of insects than under-fully lit street lights.
The two sites each have ten street lights that were sampled during 16 nights of dimmed light and during 16 nights of full light. During these total of 32 nights, we sampled insects at each street-light pole (10 per site, 20 in total). For A and also at B, we also have temperature as a variable (nightly average, max, min) - a measure per site not per light pole.
I would like to fit a model: number of  insects ~ dimming, temp_mean, temp_min, temp_max, site
I  include "site" as categorical variable and would like to use the 32 nights (16 dimmed, 16 fully lit) as 32 unique samples (each night as one sample). I am not sure whether this is "statistically sound". If not, does anybody have an idea on how to evaluate such data in order to prove the hypothesis? Would you go for t-test, PCA or something else?
streetlight <- data.frame(site=c(rep("A",16),rep("B",16),rep("A",16),rep("B",16)),
                          dimming=c(rep(1,32),rep(0,32)),temp_mean=runif(32, min=20, max=24),t
                          emp_min=runif(32, min=12, max=19),temp_max=runif(32, min=25, max=28),
                          insects=round(runif(32, min=5, max=20),0))

The hypothesis can be examined by t-test. You can group the data (number of catches) into  dimmed light, and fully lit street lights and proceed with t-test to assess whether dimmed light results  in less en-catchment or higher en-catchment or equivalent encatchment.  Alternatively, a linear model can be applied subject to certain conditions  here. PCA is not  appropriate here.

Please avoid PCA. Go with logistic regression and work out the odds ratio. This is far and away the easiest most powerful method.
Oh and you're likely to see a negative binomial distribution, very common for insect distributions. Your next step would be to get some DNA out of them