text:
Effect of mutation in DNA sequence on transcription factor binding sites
How much does a single mutation/alteration of a nucleotide affect the presence of a transcription factor binding site (TFBS)?
I am from computer science background(Obviously).
I want to make a general assumption about the number of mutated bases in a TFBS to minimize the binding affinity?

I haven't worked on TFBSs but I really, really doubt you will be able to do this. In mutations, like in real estate, it's location, location, location. It isn't a question of how many mutations but what mutations and where. TFBSs are usually described using position weight matrices, precisely because that's how you can show the importance of different changes in different positions.
Here are a couple of random Google image search results for "transcription factor binding sites logo":

(This was taken from here and you might want to read through it.)

(source: Front. Genet., 23 February 2016 | https://doi.org/10.3389/fgene.2016.00024)
The main take home messages from these images are:

Different positions have different importance and potential to disrupt binding. In the first image, I would expect that changing the T at position -3 to anything else would be pretty bad. However, changing an A to a T at -8 would likely be neutral.
Different TFs have very different binding sites and very different binding site lengths.  Look at panel B of the second image. That one appears to be a very conservative site and likely any change would be deleterious. Conversely, changing any of the nucleotides of the 17bp spacer shown in panel A would likely not be important as long as you don't change the length.

Please note that I have not done more than glance at the sources I took these from, so my inferences may well be wrong, I only use them to illustrate the general point: you cannot really put a number to this. You will always need to look at each specific change and each specific position. In some cases, a single change can disrupt binding while in others multiple changes might have no effect at all.

You might want to look into TFBS databases like JASPAR or TRANSFAC.

Every motif has a score. The score is calculated from the number of times each base was observed in that motif in the genome (I'm not sure how exactly these are calculated). Higher numbers indicate that this motif was found more frequently, which we take as a proxy for it providing better binding.
The Ensembl VEP can give you the score change for a mutated motif in the output, so this may be the easiest way to analyse this.

The motifbreakR package for R I believe does something like what you want to do;
https://www.rdocumentation.org/packages/motifbreakR/versions/1.2.2/topics/motifbreakR
If you're looking for a general rule as to what kind of mutations/where in the motif mutations will significantly affect motif binding, I suggest generating a lot of sequences with simulated SNP's (single nucleotide polymorphisms) and seeing how these score against the particular motif. Then, you will have a large(ish) database of sequences with one base pair difference from the motif and their corresponding score. This will allow you to judge what kinds of mutations affect the score the most.
Furthermore, you could somehow make a distribution of scores versus mutations and generate an empirical p-value to determine statistical significance (for instance, such and such mutation caused a score that was in the lowest 5% of scores for all mutations). You would know better than me how to do this, since your background is in CS and my PhD subject is genetics.
To make an assumption about the number of bases that need to be mutated, do this for one, two, three, etc. simulated SNP's.


text:
What score value is considered acceptable when using Vina Autodock tool for molecular docking simulations?
What score value is considered acceptable when using Vina Autodock tool for molecular docking simulations?

The answer is "Depends". But the best thing is always to look at what you ended up with and do some triaging before buying/synthesising anything. By triaging I mean filtering by its properties (.e.g. Lipinski, reactivity, synthetic feasibility etc. —discussed in this Bioinfo SE question).
Score
Regardless of docking software used, a forcefield free-energy score is a result of the scoring function used and of the properties of the molecule —including VENA which uses a hybrid score function.
This score is not, but generally is set up to try to match, the actual phys-chem one.
Gibbs free energy is measured in kcal/mol (or less commonly in kJ/mol). A hydrogen bond is about 2 kcal/mol, a π-π bond 1 kcal/mol and a salt bridge 3. The average energy of a collision of RT water molecules is 0.6 kcal/mol.
Therefore, if you dock multiple compounds (virtual library screen, VLS) you sort by free energy to see the strongest binders and consider the top ones.
However, a catch is that this will be biased by MW.
Also, it is useful to score the hits with a different program (say Gold or ICM-dock).
Ligand Efficiency
A variant metric that is employed when your hits vary in size and when they are over a given size (e.g. not fragments under 200 Da) is Ligand Efficiency (LE), which is the ratio of your score over the number of heavy atoms (ref).
Single compound
If you have a single compound and want to know if it binds acceptably, then this is a length-of-a-piece-of-string question as every compound has an affinity constant to a given protein, except most may be beyond their solubility. Initial hits will have µM affinity, while a lead will have nM affinity. However, a KD is not a ∆G and depends on your protein.

source: Stack Exchange

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

sort: votes


text:
Seeking explanation of the hg38 files downloaded from bowtie 2 website
I downloaded the H. sapiens, NCBI GRCh38 files from Bowtie's website. After unzipping, there are 6 files, 4 that end in set 1.ebwt, set2.ebwt, set3.ebwt, and set4.ebwt and two that end in set.rev.1.ebwt and set.rev.2.ebwt.
I am unable to find any explanation for what these files are either online.
which one do I use for bowtie2 mapping?
Thanks

Bowtie is providing you with the index files for hg38.  This are the result of the bowtie2-build indexer.  However you have downloaded the indexes for Bowtie1 not bowtie2

Bowtie 2’s .bt2 index format is different from Bowtie 1’s .ebwt
  format, and they are not compatible with each other.

Then according to the manual you pass the basename of the index for the reference genome. So this should be the name prior to the extension, .1.bt2 for Bowtie2.
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i>} -S [<sam>]
Main arguments
-x <bt2-idx>
The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final .1.bt2 / .rev.1.bt2 / etc. bowtie2 looks for the specified index first in the current directory, then in the directory specified in the BOWTIE2_INDEXES environment variable.

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


text:
Question on nanopore sequencing data process pipeline (cDNA-PCR)
I recently started doing the analysis on nanopore sequencing data. As I was searching for some help on pre-processing of the data, I found your nice setup pipeline created here: https://www.protocols.io/view/stranded-mapping-from-long-reads-n8cdhsw.
As I was trying to follow the instruction, I ran into a question on part where it uses lastal to divide the input file as strand specific output.
Part that uses the following script below:
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \   awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \   ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \   awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

In here, I was wondering what the ONT_barcodes_adapters.fa file actually looks like. I wanted to find the information in other places and could not find it elsewhere.
If you could help me with this, I would really appreciate it!

I recommend that you look at the most recent version of the protocol, which includes example adapter and barcode files:
https://www.protocols.io/view/stranded-mapping-from-long-reads-syheft6?step=2
[check in the protocol to see if there's a more recent version]
These files contains the public sequences that ONT have released for their sequencing kits (as briefly mentioned in the "Before Starting" section of the protocol). The sequences are available on the ONT community, but I have converted them into FASTA files for use in my protocol.
The adapter and barcode files were combined into a single FASTA file in the previous version of my protocol. For the PCR cDNA barcoding, it looked something like this:
>ONT_VNP
ACTTGCCTGTCGCTCTATCTTCTTTTT
>ONT_SSP
TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG
>BC01
AAGAAAGTTGTCGGTGTCTTTGTG
>BC02
TCGATTCCGTTTGTAGTCGTCTGT
>BC03
GAGTCTTGTGTCCCAGTTACCAGG
>BC04
TTCGGATTCTATCGTGTTTCCCTA
>BC05
CTTGTCCAGGGTTTGTGTAACCTT
>BC06
TTCTCGCAAAGGCAGAAAGTAGTC
>BC07
GTGTTACCGTGGGAATGAATCCTT
>BC08
TTCAGGGAACAAACCAAGTTACGT
>BC09
AACTAGGCACAGCGAGTCTTGGTT
>BC10
AAGCGTTGAAACCTTTGTCCTCTC
>BC11
GTTTCATCTATCGGAGGGAATGGA
>BC12
CAGGTAGAAAGAAGCAGAATCGGA
>RB12A
GTTGAGTTACAAAGCACCGATCAG

I separate out adapters and barcodes in the newer versions of the protocol because I was getting too many incorrect hits.


text:
Understanding the initialization of the Needleman-Wunsch algorithm
I would like to ask why when we initialize 0 as the starting point, there are 3 gray boxes that are not used as shown (highlighted in the red box) in the pic?

 Are these boxes to show that they are gaps or just to align them properly?

I would also like to ask how dynamic programming plays a part in filling up the rows and columns in the matrix. (I'm new to bioinformatics)

Welcome to the site Heeh. You can [edit] your question to add things to the question like more information. Could you clarify what do you mean by "how dynamic programming plays a part"? Do you want to know what is dynamic programming? Or perhaps why it is needed for this algorithm?

Hi I would like to know both about dynamic programming and why it is needed for this algorithm. I searched the meaning of dynamic programming but I do not understand the meaning behind recursive.

The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment
The dynamic programming is useful because you need first to calculate the score for each position, and then once you have the score, you can calculate the optimum path. So you solve your problem by first calculating what would happen if (there was a gap, there was a mismatch...) and then you find the optimum alignment.

Dynamic programming isn't strictly needed, it's just convenient to use recursion to handle the rows (depending on the implementation, you can also save memory by not storing the full matrix then).


text:
What should I use as a reference file for clumping?
I'm new to using PLINK and I have 2 basic questions about clumping. I'm computing polygenic risk scores for psychiatric diseases and will be using the following command to clump:
plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb

My first question: What type of file should I be using for the reference? Is it a map file?
My second question: How do I choose the most appropriate reference file?


text:
cDNA and alignment mapping
I am confused with RNA seq alignment. My understanding is that after Mature mRNA is isolated from the cell, it is then fragmented and using reverse transcriptase enzyme a cDNA copy is created which is further amplified. This cDNA is then mapped with the reference genome. If the cDNA does not have introns how does it map with the reference genome.
Can someone please help

tl;dr
Spliced aligners just make many small matches along the read pair or single read to determine splice junctions and output an alignment for the whole read, as long as it matches somewhere in the genome.
mRNA and translation
In the cell, after transcription, introns are spliced out of RNAs to become mature RNA. The protein-coding sequence does have sequence that is chimeric from cis-located portions of the same chromosome that do not actually exist.
              |<<exon 1>>|  intron    |<<exon 2>|
Genome: ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

mRNA  : ----ABCDEFGHIJKLYZ123456789aaaaaaaaaaa

As you said, if we only even consider small stretches of DNA, the sequence JKLYZ does not exist in the genome, only in the mRNA. So, we need a special way to map mRNA-seq reads to the genome that correctly splits at splice junctions.
How do we map mRNA-seq reads to a genome?
As a result, parts of the mRNA-seq reads will align well to one part of the genome, but only for part of the read. This requires a different alignment algorithm that considers multiple local high-score alignments for different parts of the read.
For example, see how a genomic aligner versus a splice aligner might treat the same mRNA-seq read on our example genome.
                         |<<exon 1>>|  intron    |<<exon 2>|
Genome:            ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

spliced aligner  : ------ABCDEFGHIJKL------------YZ123456789
genomic aligner  : ------ABCDEFGHIJKL

There is a possibility that the genomic aligner will not recognize that the second half of the read is an exon due to the large intron gap, and it may clip the alignment. Spliced aligners just take this into account and output several alignments for the same read regardless of gap size. The STAR algorithm, for example, uses many small matching pieces of a single read to determine where the splice junctions are. HiSat2 uses small approximate matches using a FM-index to find seeds for spliced alignment.


text:
In Uniprot how do I know which variants are mutations and which are polymorphic?
In Uniprot how do I know which variants are mutations and which are polymorphic?
I am looking at P38398 (BRCA1_HUMAN) specifically.
Many thanks.

What do you mean by "mutation"? Are you just talking about how frequent a variant is in a population?

This answer is very late but may still help others:
Have a look at the UniProt index of human polymorphisms and disease mutations at https://www.uniprot.org/docs/humsavar
Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.

Thanks for answering! No answer is ever late here, the objetive of this site to to amass a collection of useful Q&As, so the age of a question is irrelevant as long as it could still be asked by someone else. If, as your answer seems to suggest, you are affiliated with UniProt, it would be helpful if you could mention that since it would add more weight to your answers.


text:
Single-cell sequencing dataset has too many barcodes
I am analyzing a single-cell sequencing dataset from the website 10xgenomics, with 2000 cells. It is a BAM file and I am trying to obtain the individual cells per sample. I used the command
samtools view mtdnaAsorted.bam | grep CB:Z: | sed 's/.CB:Z:([ACGT]).*/\1/' | sort | uniq -c > reads_per_barcode_mtdnaA

to get a list of the individual bar codes and their reads, but there are 100,000 unique barcodes even with only 2000 cells. Does anyone know what could be the problem?
I know that there could be noise, but this is a problem for me because I'm trying to analyze coverage per cell, and if there are a lot of cells that contribute towards the total number of reads, the net coverage per cell is going to be really low

The BAM file you downloaded has all of the detected 10X barcodes instead of those that represent cells.  In a given 10X experiment the number of input barcodes vastly outnumbers the input cell count.  After the reaction some of barcoded beads are invariably sequenced because they are encapsulated with ambient RNA, dead cells etc. After processing the BAM file Cell Ranger infers which cells likely represent cells and generates gene expression matrices from these files.
You would need to download one of the output files that contains the list of passed barcodes such as the Gene / cell matrix (filtered) and then use the list of passed barocode sequences to filter your BAM file. Alternatively you could approximately select the same cells by selecting the top 2,000 barcodes with the highest numbers of reads.


text:
The slopes of gene networks
I have a gene network that I clustered into 8 parts. I calculated the degree distribution of the clusters using Cytoscape and the slope of the clusters range from $0.7$ to about $1.3$.
From literature (reference) usual scale free networks (real world networks) have a slope range of $2>γ>3$. What could be a reason for my network to have the smaller slope value?
Edit: The number of nodes of each cluster ranges from 431 to 187 and the number of edges range from 1009 to 263. The network is bipartite in nature.  The clusters are all very sparse which can be understood by the number of edges in the largest and smallest cluster. Are real world biological networks not always scale free. Can anyone point me towards publications that mention this fact. I am just confused about why my network doesn't conform to the norm but maybe it doesn't need to.

Could you add more information about your network? How many nodes does it have? How are them distributed in your 8 clusters? What kind of data are you using to create it (It seems like gene expression, but I'm not sure) ?

Yes, please [edit] and explain more about your network. 8 clusters suggests your network is likely very small. Is it also sparse, as biological networks often are? You also seem to be starting with the assumption that your network is scale free. Is it? It might not be, after all.

@terdon Aren't real world biological networks scale free in nature? Or that's what I have been led to believe. I have not been able to find any publications to the contrary yet. What more information do you think I can provide about the network. I have mentioned the number of nodes and edges.

Actual networks are rarely scale free.

@Devon, this seems to be the answer rather a comment. Thanks for providing the reference (I lost it/I didn't save it when I read it). @ The Last Word: Despite not being a scale-free you can still work with your clusters, but take note that the topology of your network is not to be trusted.

@Llopis I'm hoping that someone else wants to give a more thorough answer than my link :)

@Devon I think there was another paper about networks in biology, but in general it seems like people is using topologies analysis assuming too many things without (first) checking them.


text:
Compare two networks of the same genes between two species
I had a set of genes for honey bees. I found the homologues of those genes in the house fly and created a network of the honey bee genes and a network of the homologous fly genes using STRING.
I know that I could do a topological comparison between the networks and check if the same set of genes give the same hubs when in two different species. However, I can't find a single study where something similar has been done.
What biological inference can I draw from the topological parameters? For example, if the diameter of the fly network is smaller than the honey bee network, what inference could I draw?
UPDATE: It would be great if atleast some publications could be given where network comparisons have been done between different species.

So, you can give a try to this Cytoscape app that will allow to perform comparison between networks.
This app is not specifically designed for comparing networks across different species but it will enable you comparing any kind of networks.
These are the papers they published, which can give you some hints.
About your specific problem...I can't help you about what to look for when you're comparing networks from different species (not an expert in this specific topic) but as far I know you should look for something related to patterns occurring across different networks (i.e. species), network motifs, distribution of network centralities, conserved structures. In general, to compare networks, I would look for global properties (e.g. diameter, clustering coefficient...) while I would use local properties (e.g. degree, closeness, and network centralities in general) to compare the role of a node that you retrieved in all your networks.
Also, since you're looking for some literature, these links may help:

About conserved patterns of protein interaction in multiple species
About comparing genomic expression patterns across species
This one is slightly out of topic but can give you
some other ideas about comparing networks.


text:
Motif Discovery in Bipartite networks using Cytoscape
I am trying to find over-expressed motifs in a bipartite transcription factor gene network.My genes are not very well characterized, so finding functional motifs has been a challenge. I used the motif-discovery app on Cytoscape and calculated the total number of 3 node motifs on the network. In my network, where transcription factors regulate genes, the only kind of three node motif that is possible is a simple input module.

This is an example of a simple input module where a transcription factor regulates multiple genes
I further created random networks to my real world network and statistically (Student's t-test) proved that there are more number of single input modules in the real rather than the random network. However, I am just not feeling satisfied with just this one explanation because I would have like to draw some biological inference from this. I was thinking of also analysis the 4 node and 5 node motifs in the network, find the number of times any of my hub nodes occupy a single node of these motifs and maybe use that and prove the importance of my hub nodes (though hub nodes are worthy of more exploration just by virtue of their high degree). Has somebody done something like this previously. It is a given that hub nodes would be part of the maximum number of motifs just by virtue of their high degree and connection to multiple genes. So would this prove to be redundant?
Thank you in advance for any help that could be provided.

By motif I understand that you reffer to a pattern. But it is unclear if each TF has only 3 genes or any number of genes? How did you perform the t-test to compare your random networks with your found network? In general yes, it is said that highly connected genes are important and are highlighted for further studies.

@Llopis The TFs have multiple genes connected to each one. I just looked for 3 node motifs on the network. For doing the t-test, I generated 500 random networks similar to the one I have and (same number of nodes and edges) and I compared the total three node motifs in the random numbers against the number of 3 node motifs in my real network using the student's t-test.

I have recently develop an R package (sorry not a Cytoscape solution) to deal with this kind of situations of something grouping other elements, like pathways and genes.
The package is based on GSEABase of Bioconductor to store this relationships. Basically you should store a gene set for each TF and genes it is linked. I would recommend to use the name of the TF as the name of the gene set and then create a GeneSetCollection to be able to use the package.
The package provides with functionalities that allow you to calculate the distribution of genes per gene set (using genesPerPathway), as well as, pathwaysPerGene to know in how many other TF regulate each gene.
I am not sure of the biological inference you can take from this (I designed the package to test this) but you can also see if there is some specific preference for the TF to link to genes which are also regulated by other TF. This could be done by looking at the conditional probability (condPerPathways or sizesPerPathway).
Sorry the manual pages don't have examples and the vignettes are incomplete yet.
Also you could use the package tidygraph to calculate the centrality of the TF or of the genes you want.


text:
What tools can I use for a bacterial core/pan genome pipeline?
I want to perform a genome comparison on a group of isolates. I want to look into two broad groups of taxa and compare the accessory genome in each group. I have been using prokka (v1.12) and roary (v3.8.2) to do this but it appears the accessory_binary_genes.fa file is actually an untrue representation.
Note: gene_presence_absence.Rtab does contain all the full presence/ absence for accessory gene sets. Despite this Im still unhappy with the nomenclature of the gene groups [issue for another day]

[github issue 335] Your best to ignore the accessory_binary_genes.fa file. It is just for creating a quick and dirty tree with FastTree. The file itself is filtered to remove very common and not common variation to speedup the tree generation, hence the difference in numbers.
The top 5% and bottom 5% are excluded. It is truncated at 4000 genes.

I've been looking into alternative pipelines and a new software BPGA looks promising. Does anyone have experience with this?
I essentially want a tool which will give me the core and accessory gene sets, without the noise from partial gene hits.

You could use OrthoMCL which is an older and less easier to use implementation of the same workflow. I believe that will give you all the gene names back rather than truncating anything.

You can still use roary and use the information contained in the gene_presence_absence.csv file, which is complete and contains the gene IDs for each orthologous group

LS-BSR should be able to give you what you're looking for. See the article.
After you run the primary analysis, there is a simple documented workflow for splitting the pangenome into core and accessory, based on a user defined threshold. I'm the developer, so can help if you run into any problems.

Hi Jason, Thanks for the answer. Isnt LS-BSR meant to be less sensitive than roary? Also could you provide a link/ source for the workflow you mention. Thanks.

Also how does LS-BSR address genes broken due to contig borders?

This will obviously cause problems. The truncated gene will appear to be divergent or missing depending on where it is broken. One way that we address this problem is to map the raw data back against the pan-genome and determine if the region is truly missing or not, although this is not part of the LS-BSR package

I think these are really important insights that you should integrate to your answer by [edit]ing it. Also, do not hesitate to use hyperlinks to your scripts, but try to give a complete answer, so next time somebody will have similar pan genome problem, he will understand what to do now.

Roary also takes into account paralogs, so sometimes two core genes are split into different groups based on their neighbour genes and they end up with different nomenclature (group_*...). As suggested by Andrew Page in the github issue I would consider the gene_presence_absence.Rtab (this contains all the orthologous genes) and remove rows corresponding to vectors only containing 1s (core genes). In this way you will have a matrix of 1 and 0 corresponding to presence/absence of a particular accessory gene in your isolates.

Hi Sergio, yes this correct but my main concern is with gene names. If a gene is broken/ split between contigs it could end up labeled as e.g. Gene_1 and Gene_2 this would then be a mis representation of the genetic content. its partially an issue with prokka too Id imagine.

You are right, contigs border can lead to a false representation of the gene content.

Apparently there could be a new software to address these issue coming soon https://github.com/tseemann/prokka/issues/244

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


text:
How to split multifasta based on partial fasta header
I have download a batch of refseq fasta files and want to split them based on strain. this is complicated by accession numbers, project IDs etc. I downloaded the genomes via eDirect from NCBI. I know I can split the output but I want them together for some preliminary parsing (remove ambiguous and plasmid sequences etc.)
See example below.
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC

I want two files:
Clostridium botulinum BKT028387.fas and Peptoanaerobacter stomatis OBRC8.fas
My actual dataset is much larger (~70,000 genomes) and not too sure about the headers in each case.

Im thinking the following could work but Im not sure how to execute it:
Using a grep command to extract the lines but split the output.
If I download the summary list from refseq and use the names as the query example:
grep -Fwf refseq.list sequences.fa

Summary file

Peptoanaerobacter stomatis strain OBRC8, whole genome shotgun sequencing project 7,810,226 rc other DNA  NZ_AJXR00000000.1 GI:938474759

Clostridium botulinum BKT028387, whole genome shotgun sequencing project 4,134,593 rc other DNA  NZ_NOJY00000000.1 GI:1231817869

But I could edit it down to
Peptoanaerobacter stomatis strain OBRC8
Clostridium botulinum BKT028387

Am I one the right track. How do I finish this?

What operating system are you using? I asked you in your previous question as well but you didn't reply. The tools available to you and their behavior change depending on your OS.

Unix? Do you mean an actual Unix or macOS or Linux? Please tell us your actual OS, a unix-like system can be assumed if you're doing bioinformatics usually.

Yes Its is a linux system .

Great, then you should have the GNU tools. Good. Have a look at the updated answer. What does the summary list look like?

Ah! That summary file can simplify things enormously! Can you link to the entire file? Does it have all the strains and does it have each strain in exactly the same format as it appears in the fasta file header? I am afraid it won't since I've already found inconsistencies. It might be simpler to start from the beginning and download each strain separately.

For the files you have terdon, use the link and use 'send to: > complete record > file > summary > create file` https://www.ncbi.nlm.nih.gov/nuccore/?term=peptostreptococcaceae Im pretty confident they match any left behind I can curate manually.

OK, no, that summary is useless. It is unstructured and non-standard. For example, it has Peptostreptococcaceae bacterium and Uncultured Peptostreptococcaceae bacterium. I'm afraid this isn't the right way of doing what you need. It would probably be best to use a more sophisticated approach and save each sequence to the right file while downloading by also downloading the meta data for the sequence entry.

The SeqBuddy --pull_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:
$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
$: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas

This looks great. Is there a way to set this up to loop using an input list of names?

Yep, you just pass in the file name.

I had a similar problem and found a solution using awk. Since this is the first time I use awk, the syntax might not be perfect, but you should get an idea of what is possible.
awk -F' ' '{if(/^>/){sub("strain ", ""); name=$2" "$3" "$4; print > name".fa"}else{print > name".fa"}}' yourfile.fa

(more information below)
This should transform a file like this:
yourfile.fa
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis strain OBRC8 ctg120007290635 ATGC
ATGCATGC

Into two files like this:
'Clostridium botulinum BKT028387.fa'
>NZ_AESB01000237.1 Clostridium botulinum BKT028387 contig00237 ATGC
ATGCATGC
ATGCATGC
>NZ_AESB01000236.1 Clostridium botulinum BKT028387 contig00236 ATGC
ATGCATGC

'Peptoanaerobacter stomatis OBRC8.fa'
>NZ_ALNK01000035.1 Peptoanaerobacter stomatis OBRC8 ctg120007290634 ATGC
ATGCATGC
ATGCATGC
>NZ_ALNK01000031.1 Peptoanaerobacter stomatis OBRC8 ctg120007290635 ATGC
ATGCATGC

How it works

-F' ' tells awk to split the line by space (accessible with $1, $2, ...)
if(/^>/) looks for line starting with the > character, then sub("strain ", "") remove the word strain and name=$2" "$3" "$4; print > name".fa" print the line into a new file with a name based on what you asked for.
else if the line does not start with >, {print > name".fa"} appends it to a file using the current value of name.

This solution will work regardless of the number of lines between the fasta headers.
Note that I won't recommend using spaces in filenames, so I would replace name=$2" "$3" "$4 with something like name=$2"_"$3"_"$4.
I used this solution to split a large fasta file (>5 GB) containing >5k genomes with >25k scaffolds and I found it fast and rather memory efficient.

Thank you so much for the code and explanation. However, when I tried it, it gave me the following error. Any ideas?
awk: syntax error at source line 1  context is  {if(/^>/){ name=$2" "$3" "$4; print > >>>  name".fa" <<<  awk: illegal statement at source line 1 awk: syntax error at source line 1

I'm not sure about your error but you can probably put the code on several lines for easier debugging, and execute it with the awk -f option like this: awk -f myscript.awk mygenomes.fasta.

I would convert the fasta file to a format with one line per record (I have already posted the FastaToTbl and TblToFasta scripts I use for this before), and then look for a specific pattern:
FastaToTbl file.fa | grep 'BKT028387' | TblToFasta > clostridium.fa
FastaToTbl file.fa | grep 'OBRC8' | TblToFasta > peptoanaerobacter.fa

I have downloaded a few sequences using the same command you described in your previous question, and it looks like the headers are not standardized. For example, while you have many that are in the format of $species strain $strain, like these:
Clostridioides difficile strain 7496-D/ST9
Clostridioides difficile strain 5616-DH/ST42,
Peptoanaerobacter stomatis strain OBRC8
Romboutsia maritimum strain CCRI-22766
Tepidibacter mesophilus strain JCM

Others don't have the word strain and just give the strain:
Peptostreptococcus anaerobius 653-L contig00088,
Peptoclostridium acidaminophilum DSM 3953

So, as a first step, you could deal with the ones that do have the word strain as the 4th word of the header line (the 1st is the seqID which I have not included in the examples above) and the actual strain as the 5th and get those out of the way:
FastaToTbl foo.fa |
    awk '($4=="strain"){gsub("/","_"); print $0 >> $2"_"$3"_"$5".tbl"}'

That will look for lines whose 4th field is strain, then it will replace any / in the line with _ (there's a strain called CF/ST37 and that cannot be in a filename) and print into a file whose name is the 2nd, 3rd and 5th fields of the line concatenated. On my test file, this produced the following files (as you can see, some are wrong since there is a comma after the strain, but that's because they're not standardized):
$ ls *tbl
Clostridioides_difficile_5479-NonSp_ST229,.tbl  Clostridioides_difficile_6589-DH_ST42.tbl         Clostridioides_difficile_7457-NonSp_novelST,.tbl
Clostridioides_difficile_5479-NonSp_ST229.tbl   Clostridioides_difficile_6592-G_ST8,.tbl          Clostridioides_difficile_7457-NonSp_novelST.tbl
Clostridioides_difficile_5480-E_ST6,.tbl        Clostridioides_difficile_6592-G_ST8.tbl           Clostridioides_difficile_7458-N_ST10,.tbl
Clostridioides_difficile_5480-E_ST6.tbl         Clostridioides_difficile_6593-Y_ST110,.tbl        Clostridioides_difficile_7458-N_ST10.tbl
Clostridioides_difficile_5482-L_ST55,.tbl       Clostridioides_difficile_6593-Y_ST110.tbl         Clostridioides_difficile_7459-G_ST8,.tbl
Clostridioides_difficile_5482-L_ST55.tbl        Clostridioides_difficile_6595-AL_ST58,.tbl        Clostridioides_difficile_7459-G_ST8.tbl
Clostridioides_difficile_5485-Y_ST2,.tbl        Clostridioides_difficile_6595-AL_ST58.tbl         Clostridioides_difficile_7462-DH_ST42,.tbl
Clostridioides_difficile_5485-Y_ST2.tbl         Clostridioides_difficile_6596-D_ST9,.tbl          Clostridioides_difficile_7462-DH_ST42.tbl
Clostridioides_difficile_5488-DH_ST42,.tbl      Clostridioides_difficile_6596-D_ST9.tbl           Clostridioides_difficile_7463-E_ST6,.tbl
Clostridioides_difficile_5488-DH_ST42.tbl       Clostridioides_difficile_6599-NonSp_ST103,.tbl    Clostridioides_difficile_7463-E_ST6.tbl
Clostridioides_difficile_5489-M_ST15,.tbl       Clostridioides_difficile_6599-NonSp_ST103.tbl     Clostridioides_difficile_7464-Y_ST110,.tbl
Clostridioides_difficile_5489-M_ST15.tbl        Clostridioides_difficile_6600-G_ST8,.tbl          Clostridioides_difficile_7464-Y_ST110.tbl
Clostridioides_difficile_5490-CF_ST37,.tbl      Clostridioides_difficile_6600-G_ST8.tbl           Clostridioides_difficile_7466-DH_ST42,.tbl
Clostridioides_difficile_5490-CF_ST37.tbl       Clostridioides_difficile_6606-BK_ST11,.tbl        Clostridioides_difficile_7466-DH_ST42.tbl
Clostridioides_difficile_5493-DH_ST42,.tbl      Clostridioides_difficile_6606-BK_ST11.tbl         Clostridioides_difficile_7467-DH_ST42,.tbl
Clostridioides_difficile_5493-DH_ST42.tbl       Clostridioides_difficile_6607-DH_ST42,.tbl        Clostridioides_difficile_7467-DH_ST42.tbl
Clostridioides_difficile_5495-AF_ST41,.tbl      Clostridioides_difficile_6607-DH_ST42.tbl         Clostridioides_difficile_7468-NonSp_ST97,.tbl
Clostridioides_difficile_5495-AF_ST41.tbl       Clostridioides_difficile_6610-N_ST10,.tbl         Clostridioides_difficile_7468-NonSp_ST97.tbl
Clostridioides_difficile_5497-Y_ST2,.tbl        Clostridioides_difficile_6610-N_ST10.tbl          Clostridioides_difficile_7471-BI_ST1,.tbl
Clostridioides_difficile_5497-Y_ST2.tbl         Clostridioides_difficile_6612-DH_ST28,.tbl        Clostridioides_difficile_7471-BI_ST1.tbl
Clostridioides_difficile_5500-DH_ST42,.tbl      Clostridioides_difficile_6612-DH_ST28.tbl         Clostridioides_difficile_7472-Y_ST2,.tbl
Clostridioides_difficile_5500-DH_ST42.tbl       Clostridioides_difficile_6613-NonSp_ST3,.tbl      Clostridioides_difficile_7472-Y_ST2.tbl
Clostridioides_difficile_5504-L_ST55,.tbl       Clostridioides_difficile_6613-NonSp_ST3.tbl       Clostridioides_difficile_7473-DH_ST42,.tbl
Clostridioides_difficile_5504-L_ST55.tbl        Clostridioides_difficile_6614-AL_ST34,.tbl        Clostridioides_difficile_7473-DH_ST42.tbl
Clostridioides_difficile_5505-Y_ST2,.tbl        Clostridioides_difficile_6614-AL_ST34.tbl         Clostridioides_difficile_7476-Y_ST2,.tbl
Clostridioides_difficile_5505-Y_ST2.tbl         Clostridioides_difficile_6615-NonSp_ST11,.tbl     Clostridioides_difficile_7476-Y_ST2.tbl
Clostridioides_difficile_5537-D_ST9,.tbl        Clostridioides_difficile_6615-NonSp_ST11.tbl      Clostridioides_difficile_7478-NonSp_novelST,.tbl
Clostridioides_difficile_5537-D_ST9.tbl         Clostridioides_difficile_6616-NonSp_novelST,.tbl  Clostridioides_difficile_7478-NonSp_novelST.tbl
Clostridioides_difficile_5538-BK_ST11,.tbl      Clostridioides_difficile_6616-NonSp_novelST.tbl   Clostridioides_difficile_7479-NonSp_ST3,.tbl
Clostridioides_difficile_5538-BK_ST11.tbl       Clostridioides_difficile_6617-A_ST43,.tbl         Clostridioides_difficile_7479-NonSp_ST3.tbl
Clostridioides_difficile_5539-Y_ST2,.tbl        Clostridioides_difficile_6617-A_ST43.tbl          Clostridioides_difficile_7480-M_ST15,.tbl
Clostridioides_difficile_5539-Y_ST2.tbl         Clostridioides_difficile_6618-NonSp_ST53,.tbl     Clostridioides_difficile_7480-M_ST15.tbl
Clostridioides_difficile_5541-CF_ST37,.tbl      Clostridioides_difficile_6618-NonSp_ST53.tbl      Clostridioides_difficile_7481-NonSp_ST59,.tbl
Clostridioides_difficile_5541-CF_ST37.tbl       Clostridioides_difficile_6619-CF_ST37,.tbl        Clostridioides_difficile_7481-NonSp_ST59.tbl
Clostridioides_difficile_5542-Y_ST2,.tbl        Clostridioides_difficile_6619-CF_ST37.tbl         Clostridioides_difficile_7482-DH_ST42,.tbl
Clostridioides_difficile_5542-Y_ST2.tbl         Clostridioides_difficile_6621-NonSp_ST10,.tbl     Clostridioides_difficile_7482-DH_ST42.tbl
Clostridioides_difficile_5547-DH_ST42,.tbl      Clostridioides_difficile_6621-NonSp_ST10.tbl      Clostridioides_difficile_7486-A_ST43,.tbl
Clostridioides_difficile_5547-DH_ST42.tbl       Clostridioides_difficile_6622-CF_ST37,.tbl        Clostridioides_difficile_7486-A_ST43.tbl
Clostridioides_difficile_5548-Y_ST2,.tbl        Clostridioides_difficile_6622-CF_ST37.tbl         Clostridioides_difficile_7488-NonSp_ST103,.tbl
Clostridioides_difficile_5548-Y_ST2.tbl         Clostridioides_difficile_6623-NonSp_ST54,.tbl     Clostridioides_difficile_7488-NonSp_ST103.tbl
Clostridioides_difficile_5549-NonSp_ST21,.tbl   Clostridioides_difficile_6623-NonSp_ST54.tbl      Clostridioides_difficile_7491-DH_ST42,.tbl
Clostridioides_difficile_5549-NonSp_ST21.tbl    Clostridioides_difficile_6625-Y_ST2,.tbl          Clostridioides_difficile_7491-DH_ST42.tbl
Clostridioides_difficile_5550-AH_ST67,.tbl      Clostridioides_difficile_6625-Y_ST2.tbl           Clostridioides_difficile_7492-CF_ST37,.tbl
Clostridioides_difficile_5550-AH_ST67.tbl       Clostridioides_difficile_6626-Y_ST110,.tbl        Clostridioides_difficile_7492-CF_ST37.tbl
Clostridioides_difficile_5551-CF_ST37,.tbl      Clostridioides_difficile_6626-Y_ST110.tbl         Clostridioides_difficile_7494-E_ST6,.tbl
Clostridioides_difficile_5551-CF_ST37.tbl       Clostridioides_difficile_6628-AL_ST58,.tbl        Clostridioides_difficile_7494-E_ST6.tbl
Clostridioides_difficile_5552-DH_ST42,.tbl      Clostridioides_difficile_6628-AL_ST58.tbl         Clostridioides_difficile_7496-D_ST9,.tbl
Clostridioides_difficile_5552-DH_ST42.tbl       Clostridioides_difficile_6629-NonSp_novelST,.tbl  Clostridioides_difficile_7496-D_ST9.tbl
Clostridioides_difficile_5553-L_ST55,.tbl       Clostridioides_difficile_6629-NonSp_novelST.tbl   Clostridioides_difficile_7497-NonSp_novelST,.tbl
Clostridioides_difficile_5553-L_ST55.tbl        Clostridioides_difficile_6631-DH_ST42,.tbl        Clostridioides_difficile_7497-NonSp_novelST.tbl
Clostridioides_difficile_5555-DH_ST42,.tbl      Clostridioides_difficile_6631-DH_ST42.tbl         Clostridioides_difficile_7498-M_ST15,.tbl
Clostridioides_difficile_5555-DH_ST42.tbl       Clostridioides_difficile_6632-NonSp_ST3,.tbl      Clostridioides_difficile_7498-M_ST15.tbl
Clostridioides_difficile_5557-DH_ST42,.tbl      Clostridioides_difficile_6632-NonSp_ST3.tbl       Clostridioides_difficile_7499-CF_ST37,.tbl
Clostridioides_difficile_5557-DH_ST42.tbl       Clostridioides_difficile_6634-DH_ST42,.tbl        Clostridioides_difficile_7499-CF_ST37.tbl
Clostridioides_difficile_5558-AH_ST67,.tbl      Clostridioides_difficile_6634-DH_ST42.tbl         Clostridioides_difficile_7501-BI_ST1,.tbl
Clostridioides_difficile_5558-AH_ST67.tbl       Clostridioides_difficile_6635-BI_ST1,.tbl         Clostridioides_difficile_7501-BI_ST1.tbl
Clostridioides_difficile_5559-DH_ST42,.tbl      Clostridioides_difficile_6635-BI_ST1.tbl          Clostridioides_difficile_7503-Y_ST2,.tbl
Clostridioides_difficile_5559-DH_ST42.tbl       Clostridioides_difficile_6637-AL_ST58,.tbl        Clostridioides_difficile_7503-Y_ST2.tbl
Clostridioides_difficile_5573-CF_ST37,.tbl      Clostridioides_difficile_6637-AL_ST58.tbl         Clostridioides_difficile_7505-Y_ST2,.tbl
Clostridioides_difficile_5573-CF_ST37.tbl       Clostridioides_difficile_7450-DH_ST42,.tbl        Clostridioides_difficile_7505-Y_ST2.tbl
Clostridioides_difficile_5616-DH_ST42,.tbl      Clostridioides_difficile_7450-DH_ST42.tbl         Clostridioides_difficile_NCKUH-21,.tbl
Clostridioides_difficile_5616-DH_ST42.tbl       Clostridioides_difficile_7451-NonSp_ST2,.tbl      Clostridioides_difficile_RC10,.tbl
Clostridioides_difficile_5754,.tbl              Clostridioides_difficile_7451-NonSp_ST2.tbl       Clostridioides_difficile_RC10.tbl
Clostridioides_difficile_5754.tbl               Clostridioides_difficile_7452-Y_ST110,.tbl        Clostridioides_difficile_RF17,.tbl
Clostridioides_difficile_6583-AH_ST67,.tbl      Clostridioides_difficile_7452-Y_ST110.tbl         Clostridioides_difficile_RF17.tbl
Clostridioides_difficile_6583-AH_ST67.tbl       Clostridioides_difficile_7453-Y_ST2,.tbl          Peptoanaerobacter_stomatis_OBRC8,.tbl
Clostridioides_difficile_6584-AQ_ST26,.tbl      Clostridioides_difficile_7453-Y_ST2.tbl           Peptoanaerobacter_stomatis_OBRC8.tbl
Clostridioides_difficile_6584-AQ_ST26.tbl       Clostridioides_difficile_7455-A_ST43,.tbl         Romboutsia_maritimum_CCRI-22766,.tbl
Clostridioides_difficile_6588-G_ST8,.tbl        Clostridioides_difficile_7455-A_ST43.tbl          Romboutsia_maritimum_CCRI-22766.tbl
Clostridioides_difficile_6588-G_ST8.tbl         Clostridioides_difficile_7456-DH_ST42,.tbl        Tepidibacter_mesophilus_JCM.tbl
Clostridioides_difficile_6589-DH_ST42,.tbl      Clostridioides_difficile_7456-DH_ST42.tbl

Finally, you can convert these back to fasta with:
for f in *tbl; do TblToFasta $f > ${f/.tbl/.fa} && rm $f; done

This will only get those sequences whose naming scheme is as described above. But that's a start. You can then see what else you have and extract the rest manually or by finding a similar approach that matches those.

@AudileF we can't know that unless we know the details of your sequences. it all depends on how the strains are defined. I assume this is about the sequences you downloaded using edirect from your other question but that is the type of thing you should mention in your question. It's only luck that I happen to have read the other one. I am trying to give you a more general solution but the names of the few sequences I downloaded are not standard so it will not be easy.

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


text:
How do I get the gene order from genbank file using gene names?
I have a list of gene names which I know are present in a genbank file. I want to order these genes based on the order of their presence in the genbank. This should also account for whether they are on the positive or negative strand. So for example I have four genes; genes_1? are on the positive strand and genes_2? are negative
>nano Gene_names.txt
gene_1A  #assume this is the first gene for positive strand
gene_1B
gene_2A #assume this is first gene for negative strand
gene_2B

>nano Desired_out.txt
gene_1A
gene_2B
gene_1B
gene_2A

Edits
The genbank would be the standard flat file as linked below by terdon.
The gene names are from a text list file and I have edited them now to reflect this.Im not exactly sure how to approach the strandness issue but I would imagine the complement would be helpful. I am unsure as to the strandness of the genes as is.
A working example of my work is to get the order of SpoI-V and other sporulation genes from C. difficile 630
e.g. sigE, sigG, spoIIE, sleC

Please [edit] your question and clarify. I assume you mean a genbank flat file, like this one, right? If so, why does your example have a fasta-like header? Should we be looking for the "complement" field to get the strand? And should there always be one +, then one -, then one + and so on? What if there are more genes on one strand than on the other? Can you add an actual example of your data so we understand what you need?

I can't quite grok your sorting scheme, but you can use genbankr to read your genbank file into a genomic ranges object in R, and filter/sort it however you like:
source("https://bioconductor.org/biocLite.R")
biocLite("genbankr")
library(genbankr)
gb = readGenBank("sequence.gb")
subset(cds(gb), gene %in% c('sigE', 'sigG', 'spoIIE', 'sleC', 'REV7'))

output:
GRanges object with 1 range and 10 metadata columns:
      seqnames       ranges strand |        type        gene        note
         <Rle>    <IRanges>  <Rle> | <character> <character> <character>
  [1]       IX [3300, 4037]      - |         CDS        REV7        <NA>
      codon_start     product  protein_id   translation     loctype
        <numeric> <character> <character> <AAStringSet> <character>
  [1]           1       Rev7p  AAA98667.1     MNR...SLF      normal
          gene_id transcript_id
      <character>   <character>
  [1]        REV7        REV7.1
  -------
  seqinfo: 1 sequence from U49845.1 genome


text:
NCBI - edirect suite to download all genome sequences associated with a query - troubleshooting
Not sure if this is allowed but I cannot think of a better place to ask. I am attempting to download all genomic sequences from refseq associated with the query Peptostreptococcaceae. I have been been using the following command (based on some reading of the manual):
esearch -db nuccore -query "Peptostreptococcaceae" | efilter -molecule genomic -source refseq | efetch -format fasta

But I keep getting the following error
Unknown option: molecule
Unknown option: source
Must supply -query or -days or -mindate and -maxdate arguments on command
line Db value not found in fetch input No do_post output returned from
'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=nuccore&term=Peptostreptococcaceae&retmax=0&usehistory=y&edirect=3.80&tool=edirect&email=lgunn@sonic15' Result of do_post http request is .............
(a lot of text which I assume is .html or .xml related)

What operating system are you using? I just installed edirect using the instructions you linked to, then ran the exact command you show and it seems to be working fine. Could you be using an older version of edirect that doesn't support these options?

Quite possible. I tried looking up versions but could not find the info. Ill remove it and reinstall it. Fingers crossed.

That was it @terdon. Thanks. Theres no mention of updates or versions so I assumed the one I had would work fine.

Great! Could you please post a quick answer (and then accept it when the system lets you—you need to wait a bit) explaining this? That way, the question can be marked as solved.

The issue here was I was using an older version of eDirect. There is no command to check version or updates. Your best bet is if you experience problems uninstall and reinstall.

If you want to find version info you can either :
Ask for another edirect tool version (it will be the same version for edirect) :
efetch -version

6.80

Check the info contained in the edirect.pl script :
grep -i 'version = ' edirect.pl

$version = "6.80";


text:
How do I generate a color-coded tanglegram?
I want to compare two phylogenies and colour the association lines based on some metadata I have. I have been using ape cophyloplot but I have not had any success in getting the lines to colour accurately according to my data (see previous question).
Note that in my actual work flow I define the colour scheme using a palette to control the colour outcome.
I want a means to make a tangle using phylogenies which I can format. Preferably in R. I like to get an output like this:

Thanks everyone for the suggestions @b.nota answer is useful to colour according to clade groups but does not address my actual question on using metadata to colour the lines. For this answer see below:
dendextend states for color_lines
a vector of colors for the lines connected the labels. If the colors are shorter than the number of labels, they are recycled (and a warning is issued). The colors in the vector are applied on the lines from the bottom up.
# Extract labels from dendrogram on the left
labels <- dendA %>% set("labels_to_char") %>% labels

#Using a metadata table with colours create a vector of colours
labels <- as.data.frame(labels)
labels2 <- merge(labels, metadata, by.x="labels", by.y="Sample.name", sort=F)
cols <- as.character(labels2$Colours)

# Make tanglegram
tanglegram(dendA, dendC, color_lines = cols)

Side note if anyone with a phylogeny needs to convert it to a dendrogram to try this on their data see my other question.

Hi AudileF, the vector of colors for color_lines takes the order of the labels in the tree. If you want to have it follow the labels in another order, you'll need to put the colors in the initial dataframe with your metadata first, and direct color_lines to the right column of the data frame. It is hard to explain without any example.

Hi b.nota I have a workable example in my previous thread linked above, if that is of any use? Just swap rtree for rcoal, like you mentioned earlier.

I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there.
The manual mentions an argument called color_lines for the function  tanglegram():
# The `which` parameter allows us to pick the elements in the list to compare
iris_dendlist %>% dendlist(which = c(1,4)) %>% ladderize %>%
   # untangle(method = "step1side", k_seq = 3:20) %>%
   set("rank_branches") %>%
   tanglegram(common_subtrees_color_branches = TRUE)

According to the manual, the code above produces an image like this:

Try to see if you can fit your color palette in there.

To get your colors in the original order, you can make a dataframe with the labels in the right order with your metadata. You can use merge for this.
So from your example code it would be something like this:
site <- structure(list(name = structure(c(1L, 3L, 4L, 5L, 6L, 7L, 8L,9L, 10L, 2L), .Label = c("t1", "t10", "t2", "t3", "t4", "t5","t6", "t7", "t8", "t9"), class = "factor"), site = c(1L, 1L,1L, 2L, 2L, 3L, 1L, 3L, 2L, 2L)), .Names = c("name", "site"), row.names = c(NA,10L), class = "data.frame")

library(ape)
library(dendextend)

t1 <- rcoal(10)
t2 <- rcoal(10)

str(site)
# The name are factors now, so make chr first
site$name <- as.character(site$name)

# make a data.frame of your labels
labels_df <- data.frame(t1$tip.label)

#merge the 2 data.frames in the right order (hence sort=F)
colors <- merge(labels_df,site,by.x="t1.tip.label", by.y="name",all.x=T, all.y=F,sort=F)

tanglegram(t1, t2, color_lines=colors$site)

Thanks b.nota. this is something akin to a workflow I had before its good but still inaccurate to get tip labels in correct order. See my answer from that thread.  I did try this just in case and it still doesnt fix my issue. Thanks though.

As an alternative (though I realise this doesn't really answer the OPs question directly), Dendroscope and Treemap can also produce these types of plots. I don't know if they'll do coloured connections off the top of my head, but I think they can differentially colour the tree branches.
Just in case that is of use!


text:
How to extract / convert gff3 CDS sequences to multifasta
I would like to extract all the CDS from a batch of genomes. I have found a perl script from BioStars but this does not seem to work for me. I would preferably like to have a script/ method which will use the locus tag as a header.
Example
My gff files are from prokka output
ATCC0000    Prodigal:2.6    CDS 243 434 .   +   0   ID=FKKLIMLP_00001;Parent=FKKLIMLP_00001_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00001;product=hypothetical protein
ATCC0000    prokka  gene    243 434 .   +   .   ID=FKKLIMLP_00001_gene;locus_tag=FKKLIMLP_00001
ATCC0000    Prodigal:2.6    CDS 1727    2131    .   -   0   ID=FKKLIMLP_00002;Parent=FKKLIMLP_00002_gene;inference=ab initio prediction:Prodigal:2.6;locus_tag=FKKLIMLP_00002;product=hypothetical protein

If you are comfortable with python/biopython you can look into this GFF parser: https://github.com/chapmanb/bcbb/tree/master/gff

I think gffread is the best choice.

I like bedtools getfasta.  My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s.  Be aware of the -s to make sure you are pulling the correct strand.  I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations.
# bedtools getfasta

Tool:    bedtools getfasta (aka fastaFromBed)
Version: v2.26.0
Summary: Extract DNA sequences from a fasta file based on feature coordinates.

Usage:   bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:
    -fi Input FASTA file
    -bed    BED/GFF/VCF file of ranges to extract from -fi
    -name   Use the name field for the FASTA header
    -split  given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
    -tab    Write output in TAB delimited format.
        - Default is FASTA format.

    -s  Force strandedness. If the feature occupies the antisense,
        strand, the sequence will be reverse complemented.
        - By default, strand information is ignored.

    -fullHeader Use full fasta header.
        - By default, only the word before the first space or tab is used.

The gffread utility in Cufflinks package might be interesting for you. To generate a multi-fasta file with nucleotide sequences from your GFF file, then you can try:
gffread -w output_transcripts.fasta -g reference_genome.fa input_transcripts.gff

Using python and this GFF parser that mimics Biopython's SeqIO parsers:
from BCBio import GFF

# Read the gff
for seq in GFF.parse('my_file.gff'):
    # only focus on the CDSs
    for feat in filter(lambda x: x.type == 'CDS',
                       seq.features):
        # extract the locus tag
        locus_tag = feat.qualifiers.get('locus_tag',
                                        ['unspecified'])[0]
        # extract the sequence
        dna_seq = seq = str(feat.extract(seq).seq)
        # simply print the sequence in fasta format
        print('>%s\n%s' % (locus_tag, dna_seq))

AFAIK GFF3 is a version of GFF, and the page that advertises the BCBio parser indicates that it is compatible in reading/writing GFF3: http://biopython.org/wiki/GFF_Parsing

The xtractore program from the AEGeAn Toolkit was designed for this type of use case. Just set --type=CDS.
$ xtractore -h

xtractore: extract sequences corresponding to annotated features from the
           given sequence file

Usage: xtractore [options] features.gff3 sequences.fasta
  Options:
    -d|--debug            print debugging output
    -h|--help             print this help message and exit
    -i|--idfile: FILE     file containing a list of feature IDs (1 per line
                          with no spaces); if provided, only features with
                          IDs in this file will be extracted
    -o|--outfile: FILE    file to which output sequences will be written;
                          default is terminal (stdout)
    -t|--type: STRING     feature type to extract; can be used multiple
                          times to extract features of multiple types
    -v|--version          print version number and exit
    -V|--verbose          print verbose warning and error messages
    -w|--width: INT       width of each line of sequence in the Fasta
                          output; default is 80; set to 0 for no
                          formatting

It looks like you're processing prokaryotic genomes, but this program works also on eukaryotic genomes where CDSs often have introns interspersed.

You can use agat_sp_extract_sequences.pl from AGAT.
See here for all option: https://agat.readthedocs.io/en/latest/tools/agat_sp_extract_sequences.html.
Below example of use:


text:
What are the labels in my rstats distance matrix?
I have a distance matrix generated by hierfstat, thusly, with a link to the fasta file here:
library(adegenet)
library(hierfstat)
snps <- fasta2DNAbin('test.fa', chunkSize = 50)
gi <- DNAbin2genind(snps)
# manually define the populations for passing to pairwise.fst
p <- c('Botswana', 'Botswana', 'Botswana', 'Botswana', 'France', 'France', 'Vietnam', 'Vietnam', 'Uganda', 'Uganda', 'Uganda', 'Uganda', 'Vietnam', 'Vietnam', 'Laos', 'Laos', 'Laos', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam')
f <- pairwise.fst(gi, p, res.type = c('dist', 'matrix'))
as.matrix(f)

which produces
    1   2   3   4   5
1   0   0.2189008   0.225567409 0.1821518   0.259409722
2   0.2189008   0   0.130736953 0.1648034   0.191050772
3   0.2255674   0.130737    0   0.1669077   0.006396789
4   0.1821518   0.1648034   0.166907683 0   0.203931457
5   0.2594097   0.1910508   0.006396789 0.2039315   0

Is it safe to assume that 1 = Botswana, 2 = France, 3 = Vietnam, 4 = Uganda and 5 = Laos? i.e. that the distance matrix rows/columns follow the order that each population first appeared in p?
Is there a way to determine for certain, which row/column in the distance matrix corresponds to which population?

Welcome to Bioinformatics! You can change the order of your columns in the input and see how is the output, if the columns in the output follow the columns in the input, then  it would be safer. But this is from adegenet package, isn't it? Could you provide some test data to try myself? Otherwise you could check the source code of pairwise.fst

added a test file and the libraries I was using :-) i guess i could go into the source code or permute the samples to check, but seems sensible that there would be a option just to print them somehow...

Disclaimer: I don't have adagenet and hierfstat installed, I just looked at the source code. It seems like pairwise.fst (the deprecated code from the adagenet package, which you seem to be using) already returns a matrix. Set res.type="matrix" and try calling rownames(f) and colnames(f) afterwards.

Worked a treat, thanks! I'm using pairwise.fst because it is the method mentioned in this tutorial from 2015 http://adegenet.r-forge.r-project.org/files/tutorial-basics.pdf. Where did you read that it is deprecated?

I was searching for the source code. There are pairwise.fst functions in adagenet and in hierfstat. If you look at the function in the adagenet repo, you will see that it is commented out (R/fstat.R). It is, to my (limited) understanding, one of the functions that have been outsourced/exported to the hierfstat package.


text:
How to extract all pair end fastq.gz files from multiple subdirectory to a single folder
I have all file names for each paired-end fastq pair as shown below. Could you please suggest how I can extract all fastq.gz files in a single folder?
sample name:
188833984]$ ls

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc
Leaf_T1_FD_R2_L001-ds.7db8eb7e3426486db549426601b3a0bd
Leaf_T1_FD_R3_L001-ds.147177ecc03a46ccbbce3162d1185a0a

pair end file within sample name Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc
188833984]$ cd Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc

Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc]$ ls
Leaf-T1-FD-R10_S73_L001_R1_001.fastq.gz
Leaf-T1-FD-R10_S73_L001_R2_001.fastq.gz

I need all in fastq.gz format but the problem is that each sample, pair end R1 and R2 fastq reads in a single folder like this: and I need to extract all samples files in a single folder

Can you please [edit] your question and clarify. What do you mean by "extract all fastq files"? Do you just want to decompress (go from fastq.gz to fastq)? What is the final outcome you are expecting from the example you give?

I think OP wants to decompress fastq.gz files found in various subdirectories to fastq files in a single location. OP should probably decompress them where they are and create soft-links to them in their directory of choice.

@Dr.Yogesh, Do you want to copy/move all files to one folder or decompress them (as Ram RS is suggesting)?

"this file name has R1 and R1 fastq.gz": what does this mean? Did you mean to write R1 and R2?  You show a directory, Leaf_T1_FD_R10_L001-ds.3b884c360b1e4ae185408a613b90a3bc with two files, an R1 and an R2, so what is the problem? The files seem to already be in separate directories. Can you please [edit] and show us what the final result you want is?

Go to the common parent folder under which all FASTQ.GZ files reside. They could be in sub-folders under the folder, but all of them should have the folder you're in as a common parent folder. From such a folder, you should be able to run a find
find . -name "*.fastq.gz"

to get your FASTQ.GZ files. Now, there are surely better ways of doing this but here's one straightforward way.
Let's say you need to extract the FASTQs to ~/target/dir/, you'd do something like:
for f in $(find . -name "*.fastq.gz")
do
    tgt_file="~/target/dir/${f%.*}" #removed .gz from the name
    gzip -dc $f >$tgt_file #pigz -dc would be better if you have pigz
done

EDIT: OP has clarified that they wish to copy the files as such, not decompress them.
There are two ways to achieve this goal - you can either make physical copies of the files (thereby doubling the storage space you're currently using for them), or you can create symbolic links to the files (which are like shortcuts pointing to the actual files). I'll show you both ways here:
Let's assume the files are all under /fastq/files
find /fastq/files/ -name "*.fastq.gz" -exec cp {} ~/target/dir/ \; #This creates physical copies

cd ~/target/dir
find /fastq/files/ -name "*.fastq.gz" | xargs -I v_f ln -s v_f #These lines create symbolic links that point to the original files. Do not delete/move/rename the original files or these will become invalid.

Dear Ram, I want all subdirectory fastq.gz files into a single folder but files should be in fastq.gz format only, I don't want to decompress the files.

I've updated my answer. Use symbolic links instead of copying the files. As long as you don't move the original files or rename them (or any of their parent directories), your links will work fine and take much less space than if you copied the actual files.

@Dr.Yogesh Please accept the answer if it worked for you.


text:
Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements
I used the below code to do CAP analysis to link microbial data with elements, but not sure how to do statistical analysis to show which element variation is responsible most for the difference in microbial composition. I can only able to see two axis in figure so how I can get complete information about the variation based on these mineral element and individually like which mineral responsible most for microbial diversity?
    # CAP ordinate
    cap_ord <- ordinate(
        physeq = phyloseq_prop,
        method = "CAP",
        distance = bray_not_na,
        formula = ~ ~ P + Mn + Cu +Zn + Rb + Mg )

        # CAP plot
    cap_plot <- plot_ordination(
      physeq =phyloseq_prop ,
      ordination = cap_ord,
        color = "SampleType",
        axes = c(1,2)
    ) +
        aes(shape = Region) +
        geom_point(aes(colour = SampleType), alpha = 0.4, size = 4) +
        geom_point(colour = "grey90", size = 1.5) +
        scale_color_manual(values = c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
       "#AD6F3B", "#673770","#D14285","#652926", "#C84248",
      "#8569D5")
        )

        # Now add the environmental variables as arrows
    arrowmat <- vegan::scores(cap_ord, display = "bp")

    # Add labels, make a data.frame
    arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)

    # Define the arrow aesthetic mapping
    arrow_map <- aes(xend = CAP1,
        yend = CAP2,
        x = 0,
        y = 0,
        shape = NULL,
        color = NULL,
        label = labels)

    label_map <- aes(x = 1.3 * CAP1,
        y = 1.3 * CAP2,
        shape = NULL,
        color = NULL,
        label = labels)

    arrowhead = arrow(length = unit(0.02, "npc"))

    pdf("cap2.pdf", height=8, width=8)
    # Make a new graphic
    cap_plot +
      geom_segment(
        mapping = arrow_map,
        size = .5,
        data = arrowdf,
        color = "gray",
        arrow = arrowhead
      ) +
      geom_text(
        mapping = label_map,
        size = 4,
        data = arrowdf,
        show.legend = FALSE
      )

      dev.off();

     cap_ord

     Call: capscale(formula = distance ~  P + Mn + Cu +Zn + Rb + Mg, data = data)

                                  Inertia Proportion Rank
Total          2.069e+01  1.000e+00
Constrained    5.240e+00  2.533e-01   10
Unconstrained  1.545e+01  7.468e-01   97
Imaginary     -1.090e-03 -5.268e-05    1
Inertia is squared Bray distance

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3   CAP4   CAP5   CAP6   CAP7   CAP8   CAP9  CAP10
2.5733 0.8789 0.5647 0.3797 0.3037 0.1873 0.1243 0.1030 0.0800 0.0453

Eigenvalues for unconstrained axes:
 MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8
4.790 1.343 1.138 0.916 0.784 0.649 0.518 0.501
(Showing 8 of 97 unconstrained eigenvalues)

Many thanks

without some information it may be hard to answer. For example, you mention CCA but in your function call to ordinate() you are actually using a different technique method="CAP". It might help to actually see the plotted output, to understand whether the issue is in data or a difficulty of the tool. It looks like you are using phyloseq, have you looked at the documentation? https://www.rdocumentation.org/packages/phyloseq/versions/1.16.2/topics/ordinate

If you are just interested in something like a hypothesis test, you might consider using the analysis of similarity (ANOSIM): https://en.wikipedia.org/wiki/Analysis_of_similarities. You don't really need to do ordination for this, it can just take the dissimilarity matrix.

Thanks, sorry for the confusion, I wanted to ask about the CAP analysis. I just got two axis in figure explaining the total 17% variability. I tried the to capture the output of cap_ord as updated in the script but it is also not showing about complete variability explained by these all mineral element

Can I say around 25% variability explained by these variables,  from above mention data;Constrained 5.240e+00 2.533e-01

having the plot would help us to better understand your question and guide you better. Also if you could expand the data returned. Does it return any measure of variability or weight for each of your variables? (This is not shown on the output you pasted, but not sure if there is something else it isn't printed or simply it doesn't return that information). From what package is the ordinate function?

Thanks, I have updated the result

I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.
The number of constrained axes is a parameter chosen as an input to the CAPscale method:

An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for $m$? One might be concerned that the method could suffer from a poor choice for $m$. If $m$ is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if $m$ is too large relative to $N$, then a misleading canonical plot could result. For example, we cannot set $m = N - 1$, for then the squared canonical correlations would all be equal to 1.0, because we have $N$ points!

So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.
I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute  p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:

As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).

This seems to be directly implemented as the anova.cca() function in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.
For how others have approached similar questions you can see this discussion. For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). Here is the CAPscale function in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so I think that you can pass phyloseq back into vegan functions and e.g. use those analyses directly on your cap_ord object. But see that link for more info.
For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see here.

Thanks, I will try to understand these all.


text:
how to sort and join two files based on First Column Id
I have got two files large files like this in Tab-delimited format, trying to merge in R, do we need to sort the file before merging?
 otu <- read.table("asv.txt", sep="\t", header=T, blank.lines.skip=F, check.names=F)
> head(otu)
                ID                   Soil.9.S42 Soil-27.S33.L001 Soil-45.S54.L001
1 d2ec9f3b77975c0f457e4b7413b217ff         34               51               84
2 3147790f0d5a78316fb9dd64f53b9473         28               32               52
3 97aecc1f35cc1f50db507ad71dd22367         33                0                7
4 bfad6370d28182cc6304844e9bec7fb6        545              237              127
5 5fa2a987221a1d9ca416148570c18086          0                0                0
6 8e60d301122d7aa359eb6b0b00f37f62        124              143              128

> taxa <- read.table("taxonomy.txt", sep="\t", header=F ,stringsAsFactors=F,quote="")
> head(taxa)
                  ID                     V2             V3                  V4
d2ec9f3b77975c0f457e4b7413b217ff Bacteria  Cyanobacteria    Oxyphotobacteria
3147790f0d5a78316fb9dd64f53b9473 Bacteria Proteobacteria Alphaproteobacteria
97aecc1f35cc1f50db507ad71dd22367 Bacteria Proteobacteria Gammaproteobacteria
bfad6370d28182cc6304844e9bec7fb6 Bacteria Proteobacteria Alphaproteobacteria
5fa2a987221a1d9ca416148570c18086 Bacteria Proteobacteria Alphaproteobacteria
8e60d301122d7aa359eb6b0b00f37f62 Bacteria Proteobacteria Deltaproteobacteria

Now having issue with sorting file but merge command work perfectly;
newotu <- otu[order(ID),]
Error in order(ID) : object 'ID' not found
> total <- merge(otu,taxa,by="ID")

Thanks for your help.

What do you mean 'it didn't work'?

Also I'm confused - you seem to be reading the files into R and then using unix commands to merge? I prefer to do this sort of things in R - this should help https://www.statmethods.net/management/merging.html

I just use R showing the file content but wanted to do using unix command or R command.

OK - it's better to try in R IMO. Did you try follow the link I posted?

Error while doing in R, updated the command and error in the post.

Do we need to sort file in R before merging?

There was some issue with the command like header True or False, now fix it and trying again, but do we need to sort the file in R before merging?

no issue now while merging but when trying to sort file it shows error;

it doesn't need to be sorted

Thank you so much for your help. this is a very easy way to merge files in R. I also checked the file and it merges perfectly.

@user438383 or bioinfonext could one of you please post an answer with the solution you found so this question can be marked as answered?

You should check left_join() method in dplyr package. In this way you do not need to sort the data frames, as the function automatically take care of this.
For example: merged_df <- left_join(otu, taxa, by = "ID").


text:
Confirmation about PERMANOVA and Variance partitioning (permutest)
Could you please help me to understand what is the difference between Permutational multivariate analysis of variance (PERMANOVA,adonis) and  Variance partitioning
(permutest).
I understood that permutest is used to confirm the PERMANOVA result that differences in PERMANOVA is not due to dispersion.
###adonis_calculation
BC <- phyloseq::distance(phyloseq_prop, "bray")
pmv =adonis(BC ~Season*Region, as(sample_data(physeqN3), "data.frame"),permutations = 5000)
pmv

Call:
adonis(formula = BC ~ Season * Region, data = as(sample_data(physeqN3),      "data.frame"), permutations = 5000)

Permutation: free
Number of permutations: 5000

Terms added sequenti
Total         99   19.8578                 1.00000 ally (first to last)

              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
Season         1    1.8586 1.85862  12.036 0.09360  2e-04 ***
Region         1    2.4218 2.42183  15.682 0.12196  2e-04 ***
Season:Region  1    0.7522 0.75222   4.871 0.03788  2e-04 ***
Residuals     96   14.8251 0.15443         0.74656

Many thanks

Fair question, but I do think it would use some references and reasoning why you think they should (not) be the same thing. Also could be better suited for cross-validated (the statistical Stack Exchange web).

Thanks, I understand that permutest can be done to validate the permanova (adonis) result, and if permutest result is not significant then this means we can be more confident that our adonis result is a real result, and not due to differences in group dispersions.

Great, I recommend to [edit] the question. Also, your comment kind of sound you know the answer already. Is you just need a confirmation, you could really make it easier for someone to answer if you provide relevant links.

This is done; thanks!

Permutational Multivariate Analysis of Variance (PERMANOVA) assumes no distribution which means it is a non-parameteric statistical test. According to the assumption of PERMANOVA test, each group must have same dispersion in the data. Before performing any statistical test, we must check the underlying assumptions. Variance partitioning (permutest) is a confirmatory test to check the assumption (similar dispersion). You can use permutest.betadisper for permutest. After that you can perform PERMANOVA test using adonis().
Interpretation of PERMANOVA test –

If p value is small (<0.05 ) we reject the null hypothesis (No difference between the groups) otherwise we accept the null hypothesis (Significant difference between the groups).

Coefficient of determination, R-squared value indicates the variation in the distances being explained by groups under the study. Higher value indicates a better fit.

You can further look into the confidence interval for each group.
Hope this confirms the answer to your question.


text:
Mendelian randomization
I would be thankful to you if you can help with how I can use MR from my case-control study to link DNA methylation (epic array) with toxic element exposure (arsenic) and health outcome (Cardiovascular diseases).
I do have data from 125 Cardiovascular diseases and 125 non-Cardiovascular diseases individuals matched with age, sex, and smoking.I do have phenotypes files in a csv format having sex, age, smoking, disease, samples id, iAs (arsenic) values in log format and then a separate files for methylation data. we also have SNP data for these samples.
Many thanks,

Please add more detail to this question. What format are your datasets in? Can you provide any examples of that format?

added other informations, hope this is ok now.

csv and log formats are not well specified. Please show some examples of inputs and expected outputs.

thanks gringer, I don't have whole data yet but thinking how to do this analysis step by step.

I understand linkage disequilibrium (LD) versus panmixia, which the technical term for Mendelian randomisation and for humans PLINK is a population package. The best outcome is a correlation between LD and cardiovascular disease. It is however an unlikely outcome.
This really is a machine learning question and you would need to understand "feature selection" and will require very high accuracy scores >0.95 - however you could triple the size of our control group quite easily (recommended). This is a separate question and I'm happy to answer this there.

Just to address the comments. Inheritance probability, if that is the intended statistics, I personally don't think it would work because a 125 sample is low and the patients would need to be in cohorts of familial generations. LD is unlikely the explanation with humans but other non-humans species that might work (selection advantage).
My personal view is this is a straight machine learning calculation, but it's a separate question. If your sample is heavily stratified across each familial generation then Hardy-Weinberg might work, but my experience of patient cohorts is that is really unusual and there could still be confounders. Supervised learning (subset of ML) is 100% the way forward.
.... possibly use GWAS, but you've multiple features/training targets. GWAS is usually when there's one phenotype, ML would account for interaction in a way GWAS could not (it couldn't account for non-sequence features, i.e. mixed data).

Thanks, I would like to know what step by step analysis I need to perform here; like first I can look for arsenic exposure and CpG association using linear regression and then I can also look for CpG and health outcome association and then for common CpG site for arsenic exposure and health outcome I need to take for identify with the genetics variants(IVs) that affect the methylation at these site. And finally not sure if this is the correct  approach and then what I need as input for MR.

Thanks @bioinfonext Upvotes and accepts are more appreciated than thanks. I've answered above. An alternative approach e.g. via ML is a separate question.

Thanks again, still I am not clear which R packages should I use for these analysis.

GWAS: statgenGWAS package; machine learning its caret. LD is about PLINK and other autonomous packages. If you want a single package for an easy answer then shoehorn GWAS

thanks, so first I need to find CpG sites associated with arsenic exposure and then I should use these CpG site for GWAS and should also perform GWAS with disease phenotypes (CVD/noCVD). and then use these GWAS for MR?

Or do you think here I should use structural equation modelling to find CpG sites associated with arsenic exposure and disease condition as I have very small size samples and you already point out that it may not useful for MR analysis.

@bioinfonext I think its a separate question. Small sample size changes things a lot because that means no ML and no HW. I've never use SEM; I know what it does. Do consider upvoting however.


text:
Q-Q plots and the genomic inflation factor (k)
I have done the limma association analysis using below model design but not sure how to  evaluate Potential systematic biases using Q-Q plots and the genomic inflation factor (k). Could you please help with this?
#model matrix
var<-model.matrix(~BMI + as.factor(Gender) + Age +CD8T +CD4T +NK + Bcell +Mono ,data=targets2)
fit<-lmFit(mval,var)
fit2<-eBayes(fit,trend=TRUE, robust=TRUE)
probe<-topTable(fit2,adjust="BH",coef=2,num=Inf)

write.table(probe, file = "limma.output.txt")
data <- read.delim("limma.output.txt")
head(data)

               logFC       AveExpr.    t         P.Value  adj.P.Val         B
cg03546163  0.41018776  -0.8013218  7.85104192  1.13E-12    8.65E-07    12.3051291
cg17944885  -0.4884859  2.31974279  -5.7319645  6.33E-08    0.02431349  5.15237109

#For Genomic inflation factor

# For p-values, calculate chi-squared statistic
chisq <- qchisq(1-data$P.Value,1)

#Calculate lambda gc (λgc)
median(chisq)/qchisq(0.5,1)

#QQ PLOT

library(gaston)
pdf("qq.plot.pdf", width = 12, height = 12)
qqplot.pvalues(data$P.Value)
dev.off()

What have you already tried? see e.g. here and here that i found by google quickly.

thank you so much @MaximilianPress, I am not sure do we need to filter limma output based on adj.P.val before running QQ plot or not? and second I am not sure how to proceed after this step.

adj.P.val computes q-values (or applies some other correction), so its output are not really p-values. You will need to plot the p-values themselves, not a transformation.

thanks@MaximilianPress, I updated code above for QQ plot and genomic inflation factor. I hope its correct?

I don't see any obvious issues at this point but I'll let other possibly more expert people comment.

You do not need to filter the limma output based on adj.P.Val before generating the Q-Q plot / getting the inflation factor. Your code looks good to me. Is there anything unusual about the Q-Q plot or inflation factor that makes you think something went wrong?

Thanks for all help.


text:
Queries regarding MetagenomeSeq code
I am trying to understand metagenomeSeq, I will be very much thankful if you could help me to understand this:
I do have three treatment condition (T1, T2, and T3) and want to see how these treatments affect the microbial diversity in rice roots. I am interested in the pairwise comparison between them.
I used below code but I am not sure whether it is right or wrong, I do also have few queries regarding below code. I initially imported whole data into phyloseq and then filter the root data and convert into metagenomeSeq object.
root_data <- phyloseq_to_metagenomeSeq(root_M_filter)

DO I need to used normalize count for further processing? what is setting?
normalizedMatrix <- MRcounts(root_data, norm = TRUE)
settings = zigControl(maxit = 1, verbose = FALSE)

Is it correct way to make model matrix?
mod = model.matrix(~root_data$Treatment)
colnames(mod) = levels(root_data$Treatment)
res = fitZig(obj = root_data, mod = mod, control = settings)
res
An object of class "fitZigResults"
Slot "fit":
An object of class "MArrayLM"
$coefficients
                                      T1           T2          T3
71f84e7910006f22684121564206e8ca -2.8129112 -0.321366092 0.06113294
54e204fb99c80764e964456dadd6a0e5 0.5643061 0.019018800 1.01351768
55cd1fc570879d645bbf7a3642e9b0a8 -3.1159264 0.072348076 0.02205983
65f5c31e12c277aec319e2096463f9d2 -0.1355134 -0.004567206 0.53119153
7bed62f0fef250fd831dcf13bf43f4fc 1.9469670 -0.596315540 -0.35124338
                             scalingFactor
71f84e7910006f22684121564206e8ca 1.62101725
54e204fb99c80764e964456dadd6a0e5 -0.24292482
55cd1fc570879d645bbf7a3642e9b0a8 1.50462644
65f5c31e12c277aec319e2096463f9d2 0.05833637
7bed62f0fef250fd831dcf13bf43f4fc -0.56247683

What is design here?
finalMod = slot(res, "fit")$design

It is correct option for contract.make, or should I also involved some
more options?
contrast.matrix = makeContrasts(T1 - T2, levels = finalMod)
fit2 = contrasts.fit(zigFit, contrast.matrix)
fit2 = eBayes(fit2)
topTable(fit2)
                                 logFC   AveExpr         t      P.Value
3147790f0d5a78316fb9dd64f53b9473 11.262064 14.091887 16.331949
3.482815e-15
d2ec9f3b77975c0f457e4b7413b217ff 9.466087 12.212476 10.606122
6.159311e-11
7b908379d1fd2e75d84fd9aaecfd3f77 6.907182 9.784971 9.163780
1.267984e-09
d5fdc1df6dfd31acce9a34b9d8a76590 5.404459 8.373729 7.035926
1.802316e-07
19ddf8e59673c884734cf6093c672da0 5.381863 8.449413 7.011910
1.912363e-07

You get warnings, not errors. The problem is you keep trying to load ph_mgs as a dataset from a package, while it is an object you create. However, this warnings do not prevent you to analyze the ph_mgs as you wish

@DrYogesh If the problem is unrelated then I would recommend deleting this post. But you could also [edit] this one to include the relevant section (if your problem is highly related to this code).

Please use code formatting not quotation formatting for the question.