text:
Problems in creating desired phylogenetic tree with ggtree
I am working on haplotype data and want to make a tree out of haplogroups using ggtree. I have following data in newick format.
(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));

I am using following R code to plot a tree from this data,
> library("ggtree")
> tree <- read.tree(text = "(B(B2b_M112),C(C_RPS4Y711,C1a_V20,C2_M217,C2b1a_F3985,C2C1_Z1338),E(E1b1a_V38,E1b1b_M215,E2_M75),G(G_M201,G1_M342,G2a1_Z6552,G2a2a_PF3147,G2a2b2b_PF3359,G2a2b1_M406,G2a2b2a1a_U1,G2a2b2a1b_L497,G2a2b2a1c_CTS342,G2a_L30,G2b_M3155,G2b_M283),H(H1_M69,H1a_M82),I(I_M170,I1_M253,I2a1_Isles,I2a1_S21825,I2a1a_M26,I2a1b3_L621,I2a2a_M223,I2a2b_L38,I2c2_Y16419),J(J_M304,J1a_PF7257,J1a2a1a2_P58,J1a2a2_PF7264,J1a3_Z1828,J2_M172,J2a_L25,J2a1_L26,J2a1_M319,J2a1_Z1846,J2a1_Z387,J2a1_Z467,J2a1_Z6046,J2a1_Z6063,J2a1_Z6065,J2a1_Z7671,J2a1_Z7700,J2a2_PF5008, J2b_M241,J2b1_M205,J2b2_M241,J2b2_Z2456),K(K_M9),L(L_M20,L1_M22,L1a_M27,L1b_M317,L1c_M357),N(N3,N1a1_M46,N1a2_L666),O(O1a1_B384,O1b1_F2320,O1b2_P49,O2 IMS_JST0213554,O2a1_F51,O_O2a2_F525,O3_M122),Q(Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),R(R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),T(T_M184,T_L131,T_PF5633));")
> p<-ggtree(tree)
> p+geom_tiplab()

And the displayed tree looks like this

I have following problems in this tree,

It gives a low quality image
It does not show the haplotype taxa names clearly
I want to arrange haplogroups i.e nodes in the desired according to the tree given below,

i.e B haplogroups with C haplotypes and with C, E named haplotypes etc.
I have tried to solve these problems but failed to find any solution. Is there a solution?

for the tiplab issue, pls refer to the FAQ.

The above given problems were solved using manual tree construction through MEGA. MEGA allows manual construction/manipulation of trees. The newick format obtained after manually constructing the desired tree using MEGA is:
((((((((((((T_M184,T_L131,T_PF5633),T),((R_M207,M734,P224,P280xM173,R1a_M417,R1a1_M459,R1a1a_M17,R1b_M343,R1b_M73,R1b1a1b_V1636,R1b1a2_V88,R1b1b_PH155,R2_M479,R2a_M124),R)),((Q F1096,Q_L27,Q_L56,Q_M24,Q_M346),Q)),(((((((O1a1_B384,O1b1_F2320),O1b2_P49),O2IMS_JST0213554),O2a1_F51),O_O2a2_F525),O3_M122),O)),(((N3,N1a1_M46),N1a2_L666),N)),(((((L_M20,L1_M22),L1a_M27),L1b_M317),L1c_M357),L)),(K,K_M9)),((((((((((((((((((((((J_M304,J1a_PF7257),J1a2a1a2_P58),J1a2a2_PF7264),J1a3_Z1828),J2_M172),J2a_L25),J2a1_L26),J2a1_M319),J2a1_Z1846),J2a1_Z387),J2a1_Z467),J2a1_Z6046),J2a1_Z6063),J2a1_Z6065),J2a1_Z7671),J2a1_Z7700),J2a2_PF5008),J2b_M241),J2b1_M205),J2b2_M241),J2b2_Z2456),J)),(((((((((((I_M170,I1_M253),I2a1_Isles),I2a1_S21825),I2a1a_M26),I2a1b3_L621),I2a2a_M223),I2a2b_L38),I2c2_Y16419),I),((H1_M69,H1a_M82),H)),((((((((((((G_M201,G1_M342),G2a1_Z6552),G2a2a_PF3147),G2a2b2b_PF3359),G2a2b1_M406),G2a2b2a1a_U1),G2a2b2a1b_L497),G2a2b2a1c_CTS342),G2a_L30),G2b_M3155),G2b_M283),G))),((E1b1a_V38,E1b1b_M215),(E2_M75,E))),(((((((C_RPS4Y711,C1a_V20),C2_M217),C2b1a_F3985),C2C1_Z1338),C),B2b_M112),B));


text:
Why aren't clustal alignments stable under deletion of some of the sequences?
I'm new to pairwise and multiple sequence alignment in general, but I thought I understood how clustal works -- k-tuple distance is a cheap metric that's 'probably approximately' monotonic in the real global pairwise-alignment score, so we can just calculate it between all C(n, 2) sequences and do a simple UPGMA. Clustalw algorithm description here
But I'm observing odd results. Suppose I have a file with 60 sequences and run clustal on it, and observe that seqA and seqB are paired terminal child nodes. If I then delete a handful (say 10) of the sequences, but keep seqA and seqB, then in the new clustal alignment, seqA and seqB may not remain paired.
This is strange to me, because if Clustal was actually doing the k-tuple distance + UPGMA thing I thought it was doing, then all terminal pairings should always be stable under deletions of other sequences.
So... what am I missing?

Precisely that the graph of sequences depends on the sequences in the input. If you change the input sequences you change the graph, so the alignment changes too. As far as I understand the k-tuple distance has changed so the alignment too. Nice first question and welcome to the site!

@Llopis should this not be the answer, rather than a comment?

@Matt Thanks for the reminder. It might be, but I hoped that anyone could offer more insight on how the UPGMA works. Anyway I will write it.

As the graph of sequences depends on the sequences in the input, if you change the input sequences (by deleting some sequences) you change the graph.
The alignment between the sequences changes, and thus the relative distance between them changes too.
The deletion of some sequences won't cause any change on other sequences only if the full space of the C(n 2) pairwise comparisons would be done. However for a large number (and long) of sequences, this would require too much computational effort and time, making such methods as the k-tuple distance competitive.

source: Stack Exchange

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

sort: votes


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.

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


text:
Interpreting bayesian phylogenetic model using Tracer?
I am doing phylodynamic analysis of virus sequences using bayesian method using BEAUti and BEAST. As a part of my analysis I am using Tracer to check how my model is working.
I can interpret some of the results like from the Tracer like :
clock.rate - rate of nucleotide substitution
treemodel.rootheight- time when divergence occured
ESS- measures sampling quality (ideally has to be greater that 200)
But there are other many terms on this analysis which I have no idea how to explain.
Please can anybody help me to find out what other aspects of the evolution can be explained using this table.
posterior   -33531.076
prior   -481.982
likelihood  -33049.094
treeModel.rootHeight    32.271
tmrca(Brazil)   3.853
tmrca(Canada)   16.644
tmrca(China )   11.752
tmrca(Columbia) 1
tmrca(Thailand) 1.093
tmrca(USA)  32.271
tmrca(allseq)   32.271
constant.popSize    43.595
CP1.kappa   18.716
CP2.kappa   17.16
CP3.kappa   22.796
frequencies1    0.224
frequencies3    0.221
frequencies4    0.26    5001
CP1.mu  0.76
CP2.mu  1.221
CP3.mu  1.019
clock.rate  1.354E-3
meanRate    1.354E-3
CP1.treeLikelihood  -9462.258
CP2.treeLikelihood  -12520.305
branchRates 0E0
coalescent  -467.376

Did you tried reading about the tools and methods you used and what is the output you get and how to interpret it? Also could you confirm it is this tool?

I did. I followed all the tutorials and documentations. I could not find something that helps me to understand how these parameters should be interpreted.

Could you provide the links. Maybe someone reading the same could understand it better

Here is the link http://beast.community/phylodynamics_of_epidemic_influenza

Remember that you can [edit] the post to improve and add more information

You should probably pick up a basic textbook for statistical phylogenetics to get more familiar with the methods you're using. I know the BEAST people have published one (probably Drummond), but I recommend Ziheng Yang's excellent Molecular Evolution: A Statistical Approach

When I have a moment I will fill in the blanks
You must have an RNA virus (nothing else will mutate that fast). The basal group and root is the USA. It estimates the time to that basal event is 32.3 (likely years).
"frequencies" frequency of ACGT where the missing one will be 0.25.
Its a simple model of constant population size.
ESS is the most important parameter to ensure convergence of the Bayesian simulation .
Likelihood is the second most important parameter. I'll edit again shortly.

I'll fill in one more blank and summarise ...
Kappa ...
This will be the transition transversion ratio in HKY. I suspect it also will be transition transversion ratio in Kimuras 3-parameter model, this is because ... the transversions are further split into two mutation groups, ie. purine to pyrimidine has 2 further combinations and you end up with 2 ratios.  that makes sense, transition/transversion (kappa 1), transition/transversion 1 (kappa 2), transition/transversion2 (kappa 3) ....

CP1.kappa 18.716
CP2.kappa 17.16
CP3.kappa 22.796

CP2 and CP3 are either side of CP1 .. which is what you would expect on my interpretation.
These are very high values by the way, essentially there are very few transversions in this data, so ... I dunno. You need to check the data, if you see lots (and lots) of A-G and T-C mutations, then thats what I believe it is saying.
Briefly posterior and prior are the fundamental method of a Bayesian calculations, the output is probability expressed as a log, so the more negative the higher the probability (the prior looks like a random tree, or a reasonable tree with really rubbish parameters, likely the former).

Summary
Essentially I am going to call it a day here. The real issue is that the output summarise the entire modern theory in phylogenetics, so its a very long response required to address all of the parameters and I honestly don't think it helps your analysis. What I am clear about is:
* You've omitted a gamma distribution (I think Tracer would summarise this) *
So my advice is:

ESS MUST be homogeneous between runs - big trouble if not
gamma distribution MUST be run

Beast can be tricky getting gamma to converge via ESS, just a word of warning.
I can't see the full robustness, because its absent... However, from the output, I think its pretty clear, that you (well me) can more or less see the model you ran.

Thanks Michael. Yes, this is Picornavirus (RNA Virus).

It turns out there is a beast website, which I think would be worth while consulting http://beast.community


text:
Difference between copy number neutral reads and active reads
In this paper, the authors talk about copy number neutral reads (as reads that equally represent parental chromosomes) and active reads (as reads from only one parent chromosome):

We reasoned that the reads mapped into a genomic window can be partitioned in two sets: one set includes reads that equally represent parental chromosomes (copy number neutral reads); and the other set contains reads from only one parent chromosome (active reads).

Nonetheless, I have some difficulties understanding these two concepts. Can somebody try to expand on what they wrote?

Hard to answer this question, because it is not clear what you don't understand about the two concepts. Do you know what reads are? Chromosomes? I mean it is hard to guess at what level you get lost here.

A drawing of a genomic segment, with some reads mapped to it highlighting which are considered neutral and which active would be appreciated. Thanks.

So your question should be: "Can anyone make a drawing of a genomic segment, with some reads mapped to it highlighting which are considered neutral and which active?"

In genetics normally one gets a copy (allele) of a gene (or actually a chromosome) from each parent. If the alleles differ between parents, for example by a SNP, one would have one copy of each allele. With (NGS) sequencing you would get an even amount of reads for both copies. However with somatic mutations, hence mutations which occur after meiosis and gamete forming, one can lose or gain parts of the chromosome (copy number variation). In this situation, some genomic regions might be deleted or amplified. So the allele of one parent can be deleted, leaving you with only the allele of the other parent.
As requested, I have made a simple drawing of a situation where a region differs between parent alleles with SNPs. After a deletion of one of the allele only reads of the other allele are left.


text:
How do I find a list of PDB structures bound to a peptide inhibitor or peptide-like inhibitor?
I looked up existing databases but ones like iPPI-DB or enter link description here either provide non-peptide inhibitors or only sequences. I need structures of proteins bound to other proteins, the latter of which is an inhibitor of the former, and the complex has a structure deposited in the PDB database. Here's an example: PDB:3U1J
In case there is no existing database in literature, it'd be great even if I can screen for such PDB entries using some code, but I can't come up with an algorithm for it. Thank you!

Welcome to the site @simple_dozer! I'm not sure if there is a code name to know if a protein is an inhibitor, but you can start by subsetting all the structures with at least two genes

The quickest way to screen the PDB database is Advanced Search from either RCSB or PDBe.
Since you need entries with 2+ proteins you can set the number of protein chains to at least 2 in the RCSB search.  Or you can select "protein/protein complex" in Assembly composition in the PDBe search.
It's less clear how to search for inhibitors. If you'd look for the "inhibitor" word in the molecule names you'd get (together with the previous criterium) 1200+ entries.


text:
How to download SNP data from specific regions, population and positions?
I am totally new in Bioinformatics and I would like to apply my knowledge in feature selection on the tag SNP problems.
To do that, I've read a lot of papers and books in order to understand the main concepts and how the algorithms to this task work.
However, is been very hard for me to understand the organization of the huge amount of genomic data available over the internet.. HapMap, 1000 Genomes, dbSNP and so on... (!) I downloaded GB of data from different sources, but I really do not know what to do with then.
--
To be more specific in my question:
First I would like to run a Genetic Algorithm to find the tag SNP based on a simple linkage disequilibrium coverage measure. The authors of this GA used the following data

Three published SNP data sets were downloaded from HapMap (http://www.hapmap.org/) (NOT AVAILABLE ANYMORE) for evaluating the prediction accuracies achievable by our method:
• ENr112. 343 SNPs genotyped in population ASW (African ancestry in Southwest USA) on chromosome 2 from position 51512208 to 52012208. The genotypes of 83 individuals were used.
• ENr113. 380 SNPs genotyped in population ASW on chromosome 4 from position 118466103 to 118966103. The genotypes of 83 individuals were used.
• ENm013. 683 SNPs genotyped in population ASW on chromosome 7 from position 89621624 to 90736048. The genotypes of 83 individuals were used.

Besides that, the authors stated that

We represent a genotype g_i of length L by a sequence over {0, 1, 2}^L (since we assume only bi-allelic SNPs). 0 and 1 stand for the minor and major homozygous types respectively, meaning that the two alleles at that SNP locus are either 0 or 1; the major and minor alleles are determined by the frequency at which each allele appears in the population. A value of 2 denotes the heterozygous type for an SNP having two different alleles.

And also

To obtain phasing of genotypes, we used the GEVALT algorithm.

I believe they obtain the aforementioned data in genotype format (something with 0, 1 and 2 values only) and then phasing it to apply in the problem. The image below is a example of this process.

So

How and where can I download these data?
I know HapMap is not available due to a security issue but there are the FTP, but the FTP folders are really complex to understand. Can someone guide me to download data from there?

I hope I made myself clear.

I have update the entire question in order to make it clearer.

Could you clarify to which FTP folders do you refer? Could you provide a link? Also the question is about obtaining the data or about the genetic algorithm used? It seems to like two questions mixed together.

Llopis the question is about obtaining SNP data. I provide details about the paper I'm trying to reproduce in order to make easier to someone help me to find the data which could be used in this specific algorithm. About the FTP folder, I unfortunately have not enough knowledge to be more specific (I actually do not know which folder I need to download).

Maybe this helps (I could after update the question):  I study the VCF format as suggested by another user here, and I think it can be useful for me. However, I do not know how to filter the data to the specific population, chromosome and position considered in the question.

Sorry but no, it doesn't help. That this data can be used for an algorithm is not relevant to find it or download it. If you don't know a FTP site, then it is better not to mention any site at all. I'll propose some changes, review them and accept them if you think I improved your main question

Wait, you already downloaded GB of data? Are this from the HapMap project as you wished or not?

I downloaded data from this link https://mathgen.stats.ox.ac.uk/impute/1000GP_Phase3.html which is used in IMPUTE2 program (I've read in some paper that the genotype data was phased using this program) but it not seems to help.
I also downloaded the VCF content from 1000 Genomes FTP using ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/  but as I said I do not know how to filter the information.

To better help you, you need to ask one question each time (to avoid confusions like this one to future readers). One step each time. So, you already have the data you want or not? The goal of this question is to help you obtaining the data you want if you already have it we can't help because you don't need it.
If you want to ask for how to filter that information or how to run an analysis you can always ask a new question, with what you have, what is your goal and what you tried.

The HapMap website was suspended due to a security issue. Archived HapMap data is available via FTP, but it's a better idea to download up-to-date International Genome data from the IGSR website (which includes the HapMap individuals as a subset).
The format you use will be dependent on the data analysis program you're using. In most cases involving SNP analysis, the VCF format is probably what would be used.

Ok, gringer. Thanks for your reply. But it seems I need to be more specific in my answer. To do that, I copy a definition used a paper I was reading.
We represent a genotype g_i of length L by a sequence over {0, 1, 2}^L (since we assume only bi-allelic SNPs). 0 and 1 stand for the minor and major homozygous types respectively. A value of 2 denotes the heterozygous type for an SNP having two different alleles.
Which type of data do you think was used in that problem?
What I know for now is that the authors considered firstly the genotyped data and after phased the corresponding haplotypes.


text:
Network comparison of single cells (from sequencing data)
What methods can we use to compare networks (PPI, gene regulatory etc) created from single-cell gene expression (or proteomics) data?
There are methods that construct networks by comparing two conditions (e.g. from (gene) co-expression, or differential gene correlation (DGCA, Differential Gene Correlation Analysis)), and infer properties based on these networks.
What I'm interested here is the comparison of networks constructed for each single cell. The networks can be constructed solely from the data, or  using known interactions (edges) with the help of a database. These comparisons could be used for understanding or classifying cells.
This will probably require a network construction method (which may be hard because we do not have replicates for single cells), some form of network similarity measure / network alignment method. This is a very hard task in general, but here we have 1-to-1 correspondence for each node between the cells (networks), which should help.
Similar discussions at

Biostars: Question: How Do You Align Or Compare Two Networks
Content-Wise?
Researchgate: How can I measure similarity between two networks?

One related paper for single-cells is Chan et al., 2017, which uses groups of cells (created by other methods, e.g. clustering) to construct networks, and manually compares them (common edges / Venn diagram).

Could you clarify if you want a tool/methods to build those networks for each single cell or a tool/methods for comparing the cells according to the network you already have (ie assuming you already have a network for each cell)? if you ask about both it might be too broad, and that each question should be on its own to get good answers.

Well, I'm interested in both. As there seem to be no such methods yet, I think we can leave the question as is and make a new one when and if this will be answered -- and need further discussions.

Could you please clarify what type of comparison you mean. Are you only interested in common/different edges? Or are you also trying to compare network topology, shortest path analyses, centrality, network partitioning (what method?), etc etc. There are many metrics you can calculate for networks. If you just want to find the shared edges, it is trivial, assuming you have the networks. Is that really all you need?

I cannot answer this: I believe what you are asking is the essence of my question; i.e. what metrics are available, and which of these "comparisons could be used for understanding or classifying cells".

You can't build a network of a single cell only with the expression of a single cell. You either need previous known interactions or pathways or you need to use several cells/samples.
If you use previous known information, you can use pathway information, otherwise you can group some cells and use something along the lines of WGCNA to find a scale-free network.
In the first case one can add the expression upon it and see what is happening, what edges are kept and which ones not...

Once you have the network for each cell you can compare the metrics for each cell. Some metrics to look at are the centrality of each gene and the betweneess.
In a more general way, you can look if it follows a scale-free topology or other type of network structures, or how many nodes are present, how many edges are present...

Thanks for clarifying that it's not possible to build a network from data of one cell. I'd prefer using some reference interactions database anyways -- why ignore that knowledge.


text:
Different results between paired t-test, paired t-test-like (empirical Bayesian), and mixed effect model
I am learning the R package limma using its user guide. I need help on understanding the difference between the paired t-test, paired t-test-like (empirical Bayesian) and the mixed effect model in the limma.
The paired t test-like example is in the user guide 9.4.1 paired samples. The mixed effect model is in the 9.7 Multi-level Experiments (I think this is mixed effect model. Correct me if I'm wrong please).
I understand that paired t-test is a special case of linear mixed-effect modeling. So I suppose that the 9.4.1 and 9.7 should return very similar result and the 9.7 should have an identical result with the paired t-test. However, it is not the case. What am I missing here?
Please see the following small simulation study,
set.seed(20180528)

y = matrix(rnorm(2*10), nrow=2)

treat = factor(rep(c("T1","T2"),5))
subject_id = rep(c("A","B","C","D","E"), each = 2)
design <- model.matrix(~0+treat)
colnames(design) <- levels(treat)
# 1 --- 9.7 Multi-level Experiments
corfit <- duplicateCorrelation(y,design,block=subject_id)
fit <- lmFit(y,design, block = subject_id, correlation = corfit$consensus)
cm <- makeContrasts(
  Diff = T2-T1,
  levels=design)
fit2 <- contrasts.fit(fit, cm)
fit2 <- eBayes(fit2)
topTable(fit2, coef="Diff")
# logFC     AveExpr          t   P.Value adj.P.Val         B
# 1 -0.5173411 -0.08265473 -1.3646790 0.1912401 0.3824802 -4.487867
# 2 -0.1795537 -0.24626262 -0.5743662 0.5737073 0.5737073 -4.785946
# 2  --- 9.4.1 empirical bayesian paired t test
design <- model.matrix(~subject_id+treat)
fit <- lmFit(y, design)
fit <- eBayes(fit)
topTable(fit, coef="treatT2")

# logFC     AveExpr          t   P.Value adj.P.Val         B
# 1 -0.5173411 -0.08265473 -1.5408913 0.1619112 0.3238224 -4.340661
# 2 -0.1795537 -0.24626262 -0.5347977 0.6073295 0.6073295 -4.928314

# 3. -----traditional paired t test
t.test(y[1,]~treat, paired = TRUE)
# t = 1.3815, df = 4, p-value = 0.2393
# being equivalent to
# ordinary.t <- fit$coef / fit$stdev.unscaled / fit$sigma
# (1 - pt(abs(ordinary.t[1,ncol(ordinary.t)]), 4))*2

You will have better luck posting this on https://support.bioconductor.org or the limma mailing list. Also not clear why would you expect "identical results" for a paired t-test with the one using the eBayes() method

Hi, @Peter, I did post it. I expect an identical result for a paired t-test with the one in 9.7, which I think is mixed effect model. I understand that eBayes() in 9.4.1 won't have an identical result with the paired t-test.

There is a great answer on bioconductor support. Don't understand why you cross post it here...

@b.nota. Simple. Because there wasn't one.


text:
How to subset a GRanges via an argument passed into a function?
Let's say I have following example GRanges:
> library(GenomicRanges)
> gr = GRanges(
+     seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
+     ranges = IRanges(101:110, end = 111:120, names = head(letters, 10)),
+     strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
+     score = 1:10,
+     GC = seq(1, 0, length=10))
> gr
GRanges object with 10 ranges and 2 metadata columns:
    seqnames     ranges strand |     score                GC
       <Rle>  <IRanges>  <Rle> | <integer>         <numeric>
  a     chr1 [101, 111]      - |         1                 1
  b     chr2 [102, 112]      + |         2 0.888888888888889
  c     chr2 [103, 113]      + |         3 0.777777777777778
  d     chr2 [104, 114]      * |         4 0.666666666666667
  e     chr1 [105, 115]      * |         5 0.555555555555556
  f     chr1 [106, 116]      + |         6 0.444444444444444
  g     chr3 [107, 117]      + |         7 0.333333333333333
  h     chr3 [108, 118]      + |         8 0.222222222222222
  i     chr3 [109, 119]      - |         9 0.111111111111111
  j     chr3 [110, 120]      - |        10                 0
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

It is straightforward to subset this GRanges as one would subset an R data.frame, e.g.
> gr$score
 [1]  1  2  3  4  5  6  7  8  9 10

However, I am writing a function where the user would pass through a string for the metadata column to subset e.g.
subsetter = function(granges, column){
    return(granges$column)
}

I try to use the function as follows:
> subsetter(gr, score)

which outputs NULL
This obvious doesn't work, as it will subset the gr by column, not score.
> subsetter(gr, score)
NULL

If this was a data.frame, I could subset as follows:
gr[[(score)]]

but this gives the error
Error in gr[[(score)]] : this S4 class is not subsettable

Are there any other ways to subset a GRanges besides using $?

subsetter = function(gr, cname) {
    return(mcols(gr)[[cname]])
}

You can then use things like subsetter(gr, "GC") and subsetter(gr, "score").

In such cases it is best to check the str of an object gr. Then we could see that meta data is just a dataframe inside S4 object:
gr@elementMetadata
# DataFrame with 10 rows and 2 columns
#        score        GC
#    <integer> <numeric>
# 1          1 1.0000000
# 2          2 0.8888889
# ...

Once we know how to access the dataframe, we can then subset the same way as we would subset a dataframe:
gr@elementMetadata[, "score" ]

I don't see the point of wrapping this into a function, but if it is needed:
subsetter <- function(gr, cname) gr@elementMetadata[, cname ]

# test the function
subsetter(gr, "score")
# [1]  1  2  3  4  5  6  7  8  9 10

Edit
Alternatively, we could use GenomicRanges::elementMetadata() function to extract dataframe then subset:
subsetter <- function(gr, cname) elementMetadata(gr)[, cname ]

While possible, I’d generally avoid accessing S4 object slots: unless explicitly documented, they should be considered an implementation detail that may change without warning. Use the public accessors instead (in this case, mcols).

@KonradRudolph good point, added alternative solution.

Ah … well for that function the documentation says not to use it. ;-)

@KonradRudolph I see, should have read the manuals first. Thank you.

In addition to Devon’s answer, you can use non-standard evaluation to make this work with R symbol names rather than string arguments:
subsetter = function (gr, cname) {
    cname = as.character(substitute(cname))
    mcols(gr)[[cname]]
}

Now using subsetter(gr, score) works analogously to gr$score, which may be desirable. But, likewise, you’ll no longer be able to use a variable:
var = 'score'
subsetter(gr, var) # fails, same as `gr$var`


text:
Expand granges object different amounts upstream vs. downstream
I am attempting to get gene regions and their immediate neighborhoods using Bioconductor GenomicRanges packages.
It is very easy for me to extend the regions around genes if I want to extend it the same distance upstream and downstream as in this example:
library(GenomicRanges)
library(Homo.sapiens)

genes = genes(Homo.sapiens)
genes_5kb_flank = genes + 5000

And if I want to extend both upstream and downstream different amounts centered on the TSS I can use promoters(genes, upstream = 5000, downstream= 1000)
But I am having trouble finding a function that is strand-aware and also chromosome length aware that allows me to extend the gene region a different amount upstream vs. downstream... For instance if I wanted a 5kb upstream flanking region and a 1kb downstream flanking region as is used for genomic region annotations on GREAT.
Anyone know of a drop in function for this?  something like promoters() but for the full gene region instead of only around the TSS?

Can't you subset on this data? I mean do upstream +5000 and downstream +5000 and then subset for those that are above 1kb downstream. Probably you can use plyranges

There doesn't appear to be a built in function for that, so here's a modified version of promoters():
expandRange = function(x, upstream=2000, downstream=1000) {
  strand_is_minus = strand(x) == "-"
  on_plus = which(!strand_is_minus)
  on_minus = which(strand_is_minus)
  start(x)[on_plus] = start(x)[on_plus] - upstream
  start(x)[on_minus] = start(x)[on_minus] - downstream
  end(x)[on_plus] = end(x)[on_plus] + downstream
  end(x)[on_minus] = end(x)[on_minus] + upstream
  x
}

You can then call expandedGenes = expandRange(genes, 5000, 1000).

Thanks Devon this is a great.  The only problem I see with this function is that it doesn't take into account the length of the chromosomes when expanding the ranges -- when I tried to do something similar using start() and end() I kept getting errors that I expanded beyond the ranges of the chromosome.  This shouldn't be hard to fix though now that I have your function as a base, so thanks for your help!


text:
Preparing data for UpSet Plot in R
As i have multiple comparison so I m trying this UpSet plot, but it needs data to be formatted in particular way, one of the ways is as a list, so now the issue is since between one comparison the no of genes coming out may not be same as the other one, so there is a difference in dimension so if i bind them as data frame then i get dimension error, so the best way to get through is make list so what i have tried so far
A <- read.csv("GMP_MONO_UP_DOWN_CDS_EPI.txt",header = TRUE)
B <- read.csv("CMP_GMP_UP_DOWN_EPI.txt",header = TRUE)
C <- read.csv("HSC_BLAST_UPDOWN_EPI.txt",header = TRUE,sep = '\t')

dim(A)
dim(B)
ONE <- A[c(1)]#############just taking out the gene column########
TWO <- B[c(1)]#############just taking out the gene column########
THREE <- B[c(1)]#############just taking out the gene column########
dim(ONE)
dim(TWO)
b1 <- apply(ONE,2,list)
b2 <- apply(TWO,2,list)

Now this is the code from their help page
upset(movies, nsets = 7, nintersects = 30, mb.ratio = c(0.5, 0.5),
      order.by = c("freq", "degree"), decreasing = c(TRUE,FALSE))

so the data "movies" class is data frame, there is not dimension issue as such as in my case if i want to intersect two comparison which may not have the same length, so how do i overcome the issue any suggestion or help would be highly appreciated. Update my question with one more data to add.

Don't you just need a binary matrix with presence/absence data for all of the genes?
library(dplyr)
combined <- full_join(data.frame(gene=ONE, GMP_MONO=1),
                      data.frame(gene=TWO, CMP=1))
combined[is.na(combined)] <- 0

Edit:
To loop over a large number of data frames, the easiest solution is to use the reduce function from the purrr package:
combined <- reduce(list(data.frame(gene=c('gene1', 'gene2', 'gene3'), set1=1),
                        data.frame(gene=c('gene3', 'gene4', 'gene5'), set2=1),
                        data.frame(gene=c('gene1', 'gene4'), set3=1)
                       ), full_join)
combined[is.na(combined)] <- 0
combined
#    gene set1 set2 set3
# 1 gene1    1    0    1
# 2 gene2    1    0    0
# 3 gene3    1    1    0
# 4 gene4    0    1    1
# 5 gene5    0    1    0

Im getting NA in the CMP column in my combined data-frame after doing full_join , im getting 0 as well but not sure why NA..

Oh yeah. I forgot about that. See updated answer

yeah now its working but i do see error when i try to do more than ,which is apparently wont happen with full_join ,lets say i do want to intersect , more than three data frame and make a binary matrix can you suggest what more to add to your code?

see updated answer

We can plot using lists as well, using fromList function, see example:
library(UpSetR)

myGeneSets <- list(
  set1 = c("gene1","gene2","gene3"),
  set2 = c("gene1","gene4","gene5", "gene7"),
  set3 = c("gene1","gene7")
  )

upset(fromList(myGeneSets))

fromList
A function to convert a list of named vectors to a data frame
  compatible with UpSetR.

fromList(myGeneSets)
#   set1 set2 set3
# 1    1    1    1
# 2    1    0    0
# 3    1    0    0
# 4    0    1    0
# 5    0    1    0
# 6    0    1    1

another way to do the stuff cool ..i will use this as well

You need to set a common number of rows, one for each gene (make sure that the gene identifiers are the same for both files). Then in the first column you can set 1 for the file A, and 0 otherwise. In the second column you can do the same.
Then you add this dataframe as the input of the upset function.
genes <- unique(c(ONE, TWO))
df <- data.frame(A = ifelse(genes %in% A, 1, 0), B = ifelse(genes %in% b, 1, 0))
upset(df, ...)

Your solution does work but it looks like what ever intersection it does it gives only one common element ..im not sure what is wrong..

Why would you expect more than one element? Could you share the data of ONE and TWO? Otherwise could you do length(intersect(ONE, TWO))?

as I have multiple comparison between lineage type , the ONE ,TWO and now which i added later THREE they are just vector containing my gene list ,so im using only gene list ..

Then you need to add a third column and genes should be made with genes <- unique(c(ONE, TWO, THREE)), could you post your data

I will post the data ,update the question with my data.


text:
Blueprint replicates
On Encode, it's pretty easy to find out replicate information about the experiments, however I am having some difficulty finding this information on BLUEPRINT DDC, at http://dcc.blueprint-epigenome.eu
Where can I find information which describes technical replicates and biological (isogenic) replicates?
Thanks in advance.

Is that information about the samples not available on GEO?

Blueprint does not have technical replicates, you will find that very few sequencing based experiments ever have technical replicates. As biological replicates: each individual is regarded as a replicate.
This makes sense for the sort of experimental design which blueprint was designed to answer questions such as "Do individuals from group A have higher or lower H3K4me3 than individuals from groups B" or "Does cell type A have higher accessibility than cell type B". In each of these cases the within group variability to be considered is encapsulated in the difference between individuals.
Overall technical variation is assumed to be similar in all samples or at least not correlated with the contrast of interest (this may well not be true in the second example design I gave) it is safe to directly measure it, as long we know its distribution.
As an aside, one of the reasons that ENCODE got away with only two reps is that their reps are effectively technical (repeated preparations from the same cell lines) and thus there is next to no biological variation (which is what our algorithms are set up to find).

Unfortunately this information isn't as clearly laid out in Blueprint as it is in ENCODE. You can view metadata for experiments from one of the many data viewers (IHEC Data Portal, the Blueprint DCC you mentioned, DeepBlue, or the EGA dataset itself), but these tend to not contain technical/biological replicate information like ENCODE does.
You often have to dive into the data from the original publication to get that information


text:
How to use variables while calling pybedtools
I am trying to use bedtools and fetch sequences from a whole genome fasta file inside the script I get region coordinates as variables.
For example:
chr_ = "chr1"
start = 3000
end = 3402

I am not sure how to wrap my function, because the example provided by BedTool.sequence only uses coordinates in a string form.
Can you please help me to have something like this?
a = pybedtools.BedTool(chr_, start, end, from_variables=True)

You have to take your variables and make them into a single string variable.
for example:
chr_ = "chr1"
start = 3000
end = 3402
as_str = ' '.join([chr_, str(start), str(end)])
# or
as_str2 = "{} {} {}".format(chr_, start, end)

Once you have the string built you can use the example in the bedTool docs.
a = BedTool(as_str, from_string=True)

If you like an answer please accept it so future users know it worked.


text:
HIVdb program and FASTA sequence input
I am trying to use HIVdb(https://hivdb.stanford.edu/hivdb/by-sequences/) program. I have a HIV protease which I want to test with most commonly used PIs. But when I entered it in fasta sequence format to the HIVdb, it give me error code saying

Critical: There were no Protease, Reverse Transcriptase, or Integrase genes found, refuse to process.

For example, if I download FASTA sequence of 1HPV from PDB and paste it in HIVdb input sequence window as below,
PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYDQILIEICGHKAIGTVLVGPT
PVNIIGRNLLTQIGCTLNF

I get the error code. Any suggestion on how to input sequences in HIVdb?
(Note: I am new to this subject)

I am not familiar with this database or tool, so this is just a guess. Did you try to input DNA or nucleotide sequence?

The HIVdb page states the following in the second paragraph (emphasis mine):

Sequences can be entered as plain text if just one sequence is entered. Sequences must be entered using the FASTA format if multiple sequences are entered. Sequences can be pasted in the text box or uploaded using the File Upload option. The upper limit is currently 1000 sequences containing 3000 nucleotides per sequence.

The sequence in your question is a protein, so it consists of amino acids, not nucleotides. So, if you instead give a nucleotide sequence, for example this one:
>KX692212.1 HIV-1 isolate CH47 from Switzerland pol protein (pol) and gag protein (gag) genes, partial cds
CCTCAGATCACTCTTTGGCAACGACCCCTCGTCACAATAAAGATAGGGGGGCAACTAAAGGAAGCTCTAT
TAGATACAGGAGCAGATGATACAGTATTAGAAGAAATGAGTTTGCCAGGAAGATGGAAACCAAAAATGAT
AGGGGGAATTGGAGGTTTTATCAAAGTAAGACAGTATGATCAGATACTCATAGAAATCTGTGGACATAAA
GCTATAGGTACAGTATTAGTAGGACCTACACCTGTCAACATAATTGGAAGAAATCTGTTGACTCAGATTG
GTTGCACTTTAAATTTT

The program will work as expected:

For future reference, I found the KX692212 sequence above by taking your protein sequence and running a tBLASTn analysis against HIV sequences in NCBI's nr database. That was the first hit.

Yes, It worked. I didn't notice "nucleotides" emphasis.


text:
How to filter for "CAF=0.5747,0.4253"
In my gzipped VCF file I have many variants, and they report allele frequency such as:
##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

and
CAF=0.5747,0.4253

I would like to filter out the variants >= 1% allelle frequency. How do I do that?

Did you try any program or regex? How big is your VCF file ? (perhaps you could post some lines with  what you want and what you don't want) to test it

@Llopis I want all lines that the allele frequency of the alternative allele is >= 0.01. In my example, I'd keep it because 0.4253 > 0.01.

@Llopis My VCF files are big and gzipped. Direct regression expression won't work.

I agree with @Llopis, you should post a few lines of your vcf file to give a workable example. You should probably use AWK.

How do you want to deal with multi-allelic variants (e.g. A C,T) which will have >2 CAF values? Do you want all of them to be <1% or any? And as others have said, we can't really do this without a few example lines to work with. Ideally, give us a few lines including a multiallelic one and show us your desired output.

using vcffilterjdk:
you want to keep the variant if there is no INFO/CAF.
Otherwise, get All values of CAF as a String list (because it is declared as a string in the VCF header...), stream those data, skip the first value for REF , convert each string to double, you want to keep the variant if any allele is rare (< 0.01)
java -jar dist/vcffilterjdk.jar -e 'return !variant.hasAttribute("CAF") || variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf

I have not encountered this tool before, but in my experience, if Pierre suggests something it will work as advertised.  I was going to suggest a tool that I’d read about but never got around until last week: bioawk.  It was blazingly fast.

This could be done with bcftools view. To include variants where the first element in CAF is > x:
bcftools view -i 'INFO/CAF[0] > 0.1' ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz

However, it appears that CAF is encoded as a list of strings in this input so the > operator cannot be used. In addition, you may want more control on the filtering. In this case, I would suggest using a dedicated library in your favourite language for dealing with VCF files. In python, you could use pysam along these lines:
python -c "import pysam
invcf= pysam.VariantFile('ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz')

# Write to stdout
outvcf= pysam.VariantFile('-', mode= 'w', header= invcf.header)

for line in invcf:
    if float(line.info['CAF'][0]) > 0.1:
        outvcf.write(line)
outvcf.close()
invcf.close()
"

I would stay away from parsing VCF files with regex, awk, or custom text processing scripts as even if the format is quite well defined they can be quite tricky.

The first value is for the reference allele. Presumably, the OP is interested in filtering by the variant's frequency instead. So it would have to be based on elements 1 through the end of the array.


text:
KEGG FTP vs KEGG API
I was reading the KEGG plea and I found that it doesn't forbid using the KEGG API. Then, what is in the FTP server license for personal use/academic use that it is not covered by the API?
Or I could download all the database via the API?
PS: I couldn't find how to download which genes are in which pathways but it might be I didn't expend enough time with the queries.

I found the way to download the genes and pathways thanks to this answer using KEGGREST.
It seems that the API allows to fully download the database because it allows to download info, list, find, get, conv, link, ddi.
In fact in the message of 2011-05-21 it mentions that the API includes more options to download the KGML files that previously where only possible to get via the FTP site.

If you like to get pathways and corresponding genes as simple table format in a plain text file, use GeneSCF 'prepare_database' module. This can be simple solution.
FTP vs API:
The differences I see is,

quick access to the updated content via FTP. With that I assume using FTP version you will have access to weekly database updates. Also bulk download is possible and more convenient via FTP subscription.
But KEGG REST API updates once in month or two months or a year (not a recently updated contents).

Thanks for your reply, and welcome to the site! That's the kind of information I use but not what I was asking about. But GeneSCF assumes it can access the whole KEGG database via API. Could you confirm it?

It uses KEGG REST (https://www.kegg.jp/kegg/rest/keggapi.html) to access data only for specified organism. It never stores or reterives whole KEGG database. Check out [GeneSCF publication] (https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-016-1250-z) for complete details of workflow.

I've read the relevant parts of it and it didn't answer my question.

Sorry if that did not answer your question. I thought I answered "PS: I couldn't find how to download which genes are in which pathways but it might be I didn't expend enough time with the queries." this part. Good luck.

Oh, that got already answered in the other answer. Thanks, and my apologies

So, the API data is not updated at the same time that the FTP files. That's interesting: Is there any written documentation about this to read more about it?

You click those links I provided. You will get an idea about what I mentioned. There is no direct statement but if you compare those two pages, that's the information I get.

I must miss the point, I can only confirm the weekly updated KEGG FTP (on Tuesday JST) , not the update frequency of the KEGG REST API data. Anyway, thanks for your patience (I won't bother you more).

We implemented Bio2BEL KEGG package to download the KEGG database via its RESTful API

Welcome to the site Daniel. Thanks for sharing! I was reading the ComPath article precisely now. But I was not asking if it is possible to download data, but if it is possible to download all the database via the API, which apparently it is possible.


text:
How to search a specific sequence in BAM files for 10X experiment
I have to search a specific sequence in a set of cells (> 1k cells) from a single-cell experiment done with 10X genomics.
As input file I have a single bam file, and 24 fastqs, therefore each file contains information about several cells. How can I search for a specific sequence and underpin which single cells have it?
I should be able to search for a specific sequence in a bam using samtools, etc. But the issue here is how to demultiplex the data from different cells to understand which cell provide the sequence.

The bam file has bam tags which say what reads belong to what cells.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam


text:
Generating random BED intervals given constraints
Problem is to generate a random BED interval given the following constraints:

minimum start
maximum end
fixed length
maximum number of masked bases (similar to -maxN option in faSplit)
set of intervals to avoid overlap with
stay within chromosome size constraints

I can think of computationally intensive ways to code it up from scratch, but am wondering if there is a more efficient approach before re-inventing the wheel.
seq_records = {x.name: x for x in SeqIO.parse('path/to/genome.fa', 'fasta')}

def generate_random_interval(chrom, lower, upper, length, maxrep=1.0, avoid_intervals=None):
    while True:
        start = np.random.randint(lower, upper - length)
        end = start + length
        regenerate = False

        for interval in avoid_intervals:
            if (interval.start < end) or (interval.end > start):
                regenerate = True
        if 1. * seq_records[chrom].seq[lower, upper].count('N') / length > maxrep:
            regenerate = True

        if not regenerate:
            break

    return chrom, start, end

You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6).
To help with this, you could use fetchChromSizes from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest.
Given a chromosome and a maximumEnd that falls within its bounds, uniformly sample from [0, maximumEnd-minimumStart-maxLength), adding minimumStart to determine start position and maxLength as length to get the stop position.
To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into bedmap --echo-map-size operations on repeatmasked regions (4) and bedops -n 1 on "regions-to-avoid" (5) operations to do rejection sampling.
Perhaps something like:
$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1<threshold)' \
    | cut -f2- \
    | bedops -n 1 - regionsToAvoid.bed \
    | head -${SAMPLES} \
    > qualifyingCandidates.bed

This will be relatively efficient, because once the head process at the end writes the number of SAMPLES to qualifyingCandidates.bed, it will send a SIGINT signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.
Hope this helps!


text:
Viral Metagenomics
I am analyzing viral metagenomics data (Illumina Miseq) for the first time. I have used Ray (reference below) for de novo viral genome assembly before but I haven't done metagenomics analysis before.
I know that there are some tools like Metavelvet for metagenomics. But I still need to find time to learn those tools.
For now, can I do quick and dirty analysis using following method ?

De novo assemby using Ray

Create a personal database of all the viral ref genomes using NCBI database

Blast all the contigs generated by Ray against personal database.

I assume this should give me some idea about the viral genomes present in my reads. What could be the potential flaws of this method in comparison to the tools like Metavelvet that are dedicated to metagenomics?

Ray Meta: scalable de novo metagenome assembly and profiling.
Sébastien Boisvert, Frédéric Raymond, Élénie Godzaridis, François Laviolette and Jacques Corbeil.
Genome Biology (BioMed Central Ltd).
13:R122, Published: 22 December 2012

The specifics will depend on the experimental setup and data, but a few general comments...
1) Metagenomic data is often large volume and high in microbial diversity but low and uneven in coverage; single genome assemblers may suffer particularly since abundant (or closely related) species/strains will tend to be interpreted as repeats within a genome rather than deriving from different genomes, leading to assembly errors. Checkout CAMI for some recent benchmarks of mategenomic assmeblers.
2) You might want to think about filtering your database to certain viruses, e.g. you don't expect to find RNA viruses in DNA sequencing. If the sample is an animal body site (like human gut) then consider restricting to viruses found in that organism (if it's an environmental sample this may not be possible).
3) BLAST will be slow, consider DIAMOND for a drop in faster replacement for BLASTP. Even faster approaches use only specific marker genes (like Metaphlan2 or Mocat2) or k-mer methods (like kraken or clark). Again CAMI benchmarks many of these and other metagenomic taxonomic classifiers.
Otherwise, I think possibly a missing step is filtering the non-viral sequences. You haven't said if your sequencing was enriched for viruses or depleted for bacteria/host prior to sequencing but if it isn't then the vast majority of your sequences will likely be non-viral (e.g. bacterial or host). Removing or marking these early can speed up downstream analysis and remove false positives (e.g. the human genome contains integrated viruses, bacteria contain prophages etc.). Also it is normally important to look for and remove potential contaminants like lab reagents or cloning vectors.

Your answer was really insightful. Yes, I use DNAse/RNAse to deplete host DNA and the sample was enriched for viral nucleic acids.I am trying to use DIAMOND and let you know about the results. Thanks !


text:
Statistical approach to compare the SNP genotyping data among set of individuals
So, I have the genotyping data of about 650,000 SNPs for 96 individuals. I already know the Y DNA haplogroup of these individuals, so to some extent, I have a gross understanding of their ancestry.
What would be the best way to go about doing this? From what I have been reading, maybe a PCA would be a good idea? Suppose individuals 1-10 belong to a particular population (unique ancestary) then their SNP genotypes should be significantly similar to each other than say individuals 11-20 who belong to another population group. So, in theory I should get two clusters in my PCA plot.
I have the genotype of each SNP for each individual. How would I go about creating a PCA plot or any other alternative way of analysis for this kind of questions.

Way better. Now one of your question is, "how to make a PCA plot". There are the plenty of tutorials out there, have you tried one? Yes? What what was the problem with it?

What is the question?

PCA sounds like a good start to see a rough population structure among your 96 individuals, this step is half a sane thing to do, but also it can give you a hint how many distinct populations there are in your dataset. I have not personally run PCA on SNP data, but it seems that there are quite lot of tutorials online. For instance this one looks  nice.
Once you know how many categories you can expect in your population, you can use Structure to classify individuals to the categories.

Sorry if I wasn't able to make it clear, but ya that's precisely what I want to do. I want to see if the individuals cluster in any manner whatsoever. Did you mean PCA by any chance in your answer? So exactly how would I go about doing it. I know what PCA does in theory but I am not able to understand how I can create a PCA plot with the data I have in hand.

@user2887 Oh, yes. I meant PCA, fixed that already. Note that you can also edit your question, to make it clearer or to add more details.

Hey, I edited the answer (added a link to a blogpost where they show how to do PCA on SNP data). Seems that quite lot of them out there if you google "pca genotype data".


text:
Label protein strucuture residues in Chimera
I have generated proteins structure of capsid protein of two strains of viruses. I have modeled them showing Coulombic electrostatic potential described below where positive charge is blue and negative charge is red. I have opened them in two separate windows for comparison. I have following two issues:

How can I keep two molecules on the same plane when they are opened in two windows. Is there any function that lets us know which side or which angle the molecule is facing so that I can keep both molecules at same position for comparison.

I was able to label residue on the surface. But whenever I move the molecule some of the residue start moving and some remain static. I want to show that there is difference is electrostatic potential because of difference in the residue at some particular site. I have not been able to do that because those residue keep moving haphazardly?

It would be great if someone could help me with this.

You cant do this in Chimera, but in pymol you can write a script to do this.
In terms of the view, you can use get_view on the command line. This returns a section that you can copy and paste in - which will give you the same view every time. I recommend labelling the residue you want, then finding the position you want to view it from. Then using get_view to get that view, which you can use every time (just copy it into a text file) for a reference
This does presume your different molecules are aligned, but you can do that in pymol as well to make sure that when you use the output of get_view it is consistent.

https://pymolwiki.org/index.php/Main_Page this is a very useful resource, especially for instructions on making nice figures which is more complex in pymol than chimera.


text:
How to create a SAF file for virus sequence?
I am trying to create a SAF file for the virus sequence. I can download .gff3 and .gb format for the virus gene sequence but they are not supported by featureCounts. I need to create my own SAF. I tried making one but I am confused about the GeneID and Chr. The manual says Chr should match Chromosome number in the BAM or SAM file generated by align function. But viruses don't have chromosome. Could anybody help me with this please?
Here is the link of the virus sequence which I am trying to create a SAF table for.
https://www.ncbi.nlm.nih.gov/nuccore/AY386263.1

Quoted from R documentation website
featureCounts is a general-purpose read summarization function, which assigns to the genomic features (or meta-features) the mapped reads that were generated from genomic DNA and RNA sequencing.
... Both SAM and BAM format input files are accepted.
...paraphrasing
There is an in-built annotation that uses SAF format however annot.ext will override this for a custom SAF file.

According to the featureCounts documentation: "Users may also provide a GTF/GFF format annotation via annot.ext argument. But GTF/GFF annotation should only be provided as a file, and isGTFAnnotationFile should be set to TRUE when such a annotation is provided."

Cross-post here https://www.biostars.org/p/318489/


text:
Is there any tools that can generating FBA model from genome sequences or proteome sequences?
Suppose I have all sequences of genes or proteins of certain species, are there any free tools I can use to generate a flux balance analysis model out of them?
What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

To generate flux balance, one should now the stoichiometric numbers of all the reactions that are taking place, as well as the conditions (that could affect the reactions). I am not sure we have enough knowledge, time and computationally power to do this.
At the moment I'm only aware of a paper simulating the Mycoplasma genitalium cell. But I'm not sure it was worth the effort as I haven't seen replicated or redone with other species, perhaps the model was not accurate enough to be worth it.
The tools used are here. The licenses are not a standard, but at least they seem free.

Thanks. What I was thinking is more like to annotate all the enzymes based on some annotation method. Then mapping these enzymes to known meta-metabolic networks that we know for instance the meta-metabolic network in KEGG database.

Which would be your input? A list of genes? This is trivial to do, but you don't know how many active sites are in an enzyme nor if there is any inhibitor, and you would need to know the initial concentrations to see if you reach a steady state or not, and what is your end point is also a variable you should define


text:
Determining position of side chain hydrogen in glycine residues from coordinates of backbone atoms
I have a dataset of pdb-format protein structure decoys/models generated by tertiary structure prediction programs. Many of these are missing hydrogen atoms.
Specifically, I'd like to determine the rough location of side-chain hydrogens in glycine residues, given the CA, C, N, and O coordinates. Is this possible? I don't need too much precision - I just need the general direction relative to the rest of the residue. More specifically, I want a vector originating at the CA that points roughly in the direction of the side chain.
I've looked at software for generating hydrogen atoms from pdb structures, however, I need to know which ones are glycine side chains, and was wondering if there might be a simpler way to predict hydrogen locations given coordinates of backbone atoms.

Glycine has a single hydrogen atom as its side chain:

All the six bond angles with the CA atom in the middle are about 109°
(C-CA-N, C-CA-HA3, C-CA-HA2, N-CA-HA3, N-CA-HA2 and HA3-CA-HA2 using the CCD naming convention). This defines the rough direction of the hydrogen atoms.
But which of the two hydrogens is the side chain?
If it is a left-handed protein, I'd assume that it's the hydrogen pointing toward the screen, similarly to the side chain in L-alanine:

The images are taken from WikiMedia Commons:
glycine
and L-alanine.

Thank you! So, for clarification, if I find the two vectors where bond angles N-CA-HA and C-CA-HA are about 109 degrees, and choose the one corresponding to a left-handed protein (assuming my structures are left-handed), I will likely have a reasonable estimate for the direction vector I'm looking for?

@JonathonByrd yes, I'd think so


text:
Issues installing Dendropy
I want to perform phylogenetic analysis using Dendropy. I have installed dendropy using command sudo pip install -U dendropy but when I import it in python3.5.1 it gives followinng error:
 > python3
   Python 3.5.1 |Anaconda 2.5.0 (64-bit)| (default, Dec  7 2015, 11:16:01)
   [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] on linux
   Type "help", "copyright", "credits" or "license" for more information.
 >>> import dendropy
 Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
 ImportError: No module named 'dendropy'

I again installed dendropy using the above given command as follows:
sudo pip3 install -U dendropy
Requirement already up-to-date: dendropy in /usr/local/lib/python3.5/dist-packages (4.4.0)
Requirement not upgraded as not directly required: setuptools in /usr/local/lib/python3.5/dist-packages (from dendropy) (39.2.0)

After doing this again the error remains the same and I am still not able to import dendropy, any solution to this problem?

You should never install python packages as root or with sudo.
You're already using conda, so conda install -c bioconda -c conda-forge dendropy.

And why not install python packages as root ?

You should always do as little as root as possible to minimize (A) the odds of screwing something up and (B) possible security problems. Further, later in your career you're likely to run into cases where you don't have the root password and IT will have you fired if you try to install packages as root on a central server. pip has a --user option for exactly such reasons.


text:
Use of heterozygous SNPs in cancer research: why?
When reading about allelic fraction (AF) and SNPs in cancer research, they always mention the fact that they're using heterozygous SNPs (informative SNPs). Why is this? Why can't we use homozygous SNPs?
In this case, what is the reference base (from the reference genome)?

I’m no longer working in tumour sequencing so I’m by no means an expert. But in a nutshell, the reason is that, as indicated, homozygous SNPs aren’t informative: if your allelic fraction is 100%, meaning you only observe a single nucleotide at a given position, we don’t know whether we’re dealing with a tumour related mutation. In fact, we probably don’t because the sample likely contains healthy as well as tumour cells.
By contrast, tumour-specific variants always have an AF < 100%, in other words they exhibit some degree of heterozygosity.
In fact, in healthy tissues you’d expect to see one of three cases: your variant allele (= sequence base ≠ reference) has a frequency of either 0% or 100% (homozygous), or 50% (heterozygous), ignoring rare somatic mutations. With tumour samples, by contrast, the fraction will vary depending on the penetrance of the variant in the tumour, and on the fraction of healthy to tumour cells in the sample.

I'm not sure what I'm talking about, but I guess one could re-phrase Konrad's answer in terms of information theory.
The information content of a message (a packet of bits) is inversely proportional to the frequency of those bits. For example, if you randomly extract a word (the message) from an English book and the word is "the", you don't get much information about the book because "the" is very frequent and you see it everywhere. If, instead, the word happens to be the much rarer "bioinformatics" you start having a pretty good guess what the book is about.
About SNPs, homozygote reference genotypes have very little information because most of the genomic positions in most of the people are homozygote reference. Heterozygote and homozygote non-reference are rarer and thus more informative.
So, I don't think there is anything wrong with using homozygote SNPs. They are ignored because they are so common that their information content is very small.


text:
Disk space error while aligning reads using STAR
Hi am trying to align RNA-seq reads using STAR through following commands,
STAR --runMode alignReads --runThreadN 26  --outSAMtype BAM  SortedByCoordinate Unsorted --sjdbOverhang 100  --genomeDir wheat --readFilesIn G3_cleaned_R1.fastq G3_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3
Jun 05 13:53:59 ..... started STAR run
Jun 05 13:53:59 ..... loading genome
Jun 05 14:07:12 ..... started mapping
Jun 05 18:42:47 ..... started sorting BAM

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin   does not agree with the actual size on disk: 509186583   0   45

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 648593301   0   47

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 660412146   0   46
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 470931698   0   44

Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 653155717   0   43

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Segmentation fault (core dumped)

Although there is much free space on my system as shown below:
> df -Th
  Filesystem     Type      Size  Used Avail Use% Mounted on
  devtmpfs       devtmpfs   63G     0   63G   0% /dev
  tmpfs          tmpfs      63G  136K   63G   1% /dev/shm
  tmpfs          tmpfs      63G  1.9M   63G   1% /run
  tmpfs          tmpfs      63G     0   63G   0% /sys/fs/cgroup
  /dev/sda3      ext4      3.6T  1.4T  2.2T  39% /
  /dev/sda2      ext4      9.8G   57M  9.2G   1% /boot
  tmpfs          tmpfs      13G     0   13G   0% /run/user/480
  tmpfs          tmpfs      13G   16K   13G   1% /run/user/1000

How this problem can be solved?

Are you an older version of the STAR aligner? A similar error has been discussed here: https://groups.google.com/forum/#!topic/rna-star/OTUG26gKKmo

You are using 26 threads, while you have 125GB of RAM, this means per thread you have 4.8 GB RAM. Your genome size is 13 GB (according to your previous post). An answer to your previous post suggested to reduce the number of threads. Did you try to run with lower number of threads?

@b.nota I ran the job discussed in previous post with same number of threads.

@Kasper I am using latest version of STAR

Report the bug on github and sort with samtools.

@DevonRyan ok I will post it there. Can I use --limitBAMsortRAM to solve this problem?

@AmmarSabirCheema I don't think so, since any time you need to increase that there's an error message saying so. I think this is just a bug somewhere. Just sort the unsorted output with samtools or sambamba.

STAR versions 2.6.0b and 2.6.0a are unstable,  as written Issues performing variant calling with GATK you are using version 2.6.0b. You should switch to version 2.6.0c which is stable.  Also avoid to perform multiple alignments in the same directory.


text:
Test to determine if two genes/exons share the same evolutionary histories?
In classic phylogenetic inference one is usually given various orthologue sequences of a given gene across various species. Those sequences are then multiple aligned and used to construct a phylogenetic tree.
Say, I split the multiple alignment of this gene into two parts (part A and B).
My question now is: What is the most appropriate test to determine if part A and B share the same phylogenetic history?
Having looked through literature, it appears this can be done by constructing phylogenetic trees separately for part A and B and test for incongruence using

Kishino–Hasegawa (KH test);
Shimodaira–Hasegawa (SH-test);
SOWH-test etc.

as nicely summarized in Plantet et all, and Goldman et al. etc.
However the above test rely on constructed trees beforehand which might be incorrect due to inappropriate algorithms/parameters.
I was therefore wondering if there are more appropriate test maybe without having to construct trees at all?

I don't really see how that would be possible, off the top of my head. After all, when you construct trees, you are "operating directly on the sequences" and trees are the tool for determining evolutionary history. Could you [edit] your question and explain more about why you don't want to use trees? Why don't you construct the tree from your alignment?

Well I didn't mean to say that I don't want to use trees, just if there are methods that don't need to compute them. I was just concerned that the way the tree is constructed might introduce a bias. I'd be happy to learn that SH-test and the like are appropriate approaches to answer the stated question. I just wanted to make sure I didn't overlook more appropriate approaches.

I'm not sure that you can use any of the tests you've listed on fundamentally different datasets - the likelihoods will always differ in incomparable ways. One way you could do this is, if you have a priori two different hypotheses, calculate the site-specific likelihood of the whole alignment on each of the two topologies and compare the differences in likelihood across sites. If one segment of the alignment is biased towards one resolution, the difference in likelihoods between the two topologies should show that.

Although I see in the Plantet paper the two topological-based tests, which I guess would work. The only question I have with those is how they accommodate uncertainty. It may be, for example, that in partitioning your dataset you arrive at block B containing very little phylogenetic signal for any hypothesis, and treating its point estimate as concrete is lending unjustified support to it.

Good points @NatWH. Accommodating uncertainty is indeed an issue and I suppose one has to show that there is a sufficient amount of information bearing sites for each partitions, however I couldn't find suitable numbers anywhere. I'd be curious on any good practice guidelines or maybe papers that have done it nicely since this seems of general interests.

As for your first point @NatWH, I had the following in mind: Calculating the maximum likelihood trees (tree A and B) for both partitions (part A and B) separately. I'd then calculate the likelihoods and perform a incongruence test only for part A data under both tree topologies of tree A and B. Same procedure for part B as well for completeness. So the test is in using the same dataset(s) and the likelihood should be comparable, do you agree?

Yes I agree, that shouldn't violate the assumptions of the SH test, although I am just not familiar enough with the SOWH to say anything about that. I would still advise determining the site-specific differences in log-likelihood, though. Its helpful to know the distribution of the signal.

Agreed! I just read you are going to start at the Plant Science Dept in Cambridge, where I also work. Small world! We could meet to discuss in person at some point, I'm in the Baulcombe lab. Thanks again for your input anyway!

Absolutely! I look forward to that in September.

Since there hasn't been any other answers, I'd thought to share my approach in case someone is having the same issue.
Gauging from the comments of tendon and NatWH, it seems using the SH-test valid stategy to tackle this problem.
Having done some research, I could only find one functional software package, the R-package phangorn which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:
Reading in multiple-alignment (MLN). Here: amino-acids (AA)
 library(phangorn)
 mln <- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")

Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)
breakpoint <- 100
 # partA: 1:breakpoing
 ( mln_partA <- subset(mln, select = 1:breakpoint, site.pattern = FALSE) )
 # 58 sequences with 100 character and 45 different site patterns.
 # The states are a r n d c q e g h i l k m f p s t w y v

 ## partB: breakpoint:end
 ( mln_partB <- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.

Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method
 ## Note, other tree building algorithms might be used
 ( treeNJ_partA <- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB <- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.

Using topology of previously constructed trees and compute optimal likelihood based on partA data
 ( treeNJ_partA_fit1 <- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --> -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 <- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --> -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292

Performing SH-test for for both trees
 ( shPart1_t12 <- SH.test(treeNJ_partA_fit1, treeNJ_partB_fit1, data = mln_partA, B = 1e5) )
 #      Trees      ln L Diff ln L p-value
 # [1,]     1 -1379.041     0.000 0.48461
 # [2,]     2 -1530.682   151.641 0.00000

As can be seen in this example, the tree constructed for partB doesn't fit partA data as much as the tree constructed from partA.

Thanks for sharing! But perhaps it would be clear if the output is removed, or if the code is explained before the block

Hmm, in fact I included the output to make it clearer since I like to know intermediate results. I'll exclude a few lines though and maybe someone else can share their opinion.

Well, it is only my opinion! But for me it would be clearer if you explained the code alongside breaking the block into smaller ones. As is I find hard to  know what is going on

Agreed. I broke it into chunks, hope this improved it!

Nitpicking: whatever model dist.ml uses for amino acids by default may not be the optimal model for the dataset (given a tree), unless you selected that elsewhere and did not show it in the code. Likewise I don't think the NJ tree of ML distances is guaranteed to give the ML tree. It appears that optim.pml does not perform topological searches by default. However, the assumptions of the SH-test call for the ML tree of the data. I would instead use an external program (RAxML, IQ-TREE), to find the overall ML topology and the ML topology of each subset, and run the SH-test on those.

You could also use one of the implementations of model finder (AUTO in RAxML or MFP in IQ-TREE) to find the best-fitting model of substitution.

That's true. In fact I've used RAxML for the original data but used non-optimized NJ deliberately to make a more reproducible example and keeping the code simple. In fact, this underlines my original point: Testing incongruence between partitions shouldn't depend (too much) on the tree building algorithm or/and optimizations. If it does, I'd love to read about guidelines, best practices or examples on how to do this properly.


text:
Macs2 peak calling?
I have paired end ChIP-seq data with 101 bp and 2 biological replicates for each one. I have done peak calling with macs2 but I have some questions about it.
I also faced with an warning:

WARNING @ Thu, 07 Jun 2018 17:06:05: #2 Since the d (197) calculated from paired-peaks are smaller than 2*tag length, it may be influenced by unknown sequencing problem!
  WARNING @ Thu, 07 Jun 2018 17:06:05: #2 You may need to consider one of the other alternative d(s): 197
  WARNING @ Thu, 07 Jun 2018 17:06:05: #2 You can restart the process with --nomodel --extsize XXX with your choice or an arbitrary number. Nontheless, MACS will continute computing.

I have added --nomodel --extsize 197 ; --nomodel --extsize 147 and --nomodel --extsize 202 (separately to macs2 command) and got the results without any warning? Which one is more correct?
Are broad peaks extended of narrow peaks? If I apply intersect between them I should expect find 100% overlap between narrow peaks and broad ones?
Which kind of peak (narrow/broad) is proper for H3k27ac, H3k4me1, H3k4me3, H3k27me3 study?
If there is no control group for using as background, can I use default parameters?

cross posted: https://www.biostars.org/p/319418/

Your original command without --nomodel --extsize ... is probably the most accurate. This warning stems from a time when reads were much much shorter and likely never made that much sense to begin with.
Broad peak calling in MACS2 basically works by finding a bunch of nearish narrow peaks and merging them. If you have really broad signals then use something like histoneHMM instead.
H3K27ac, H3K4me1 and H3K4me3 are narrow. H3K27me3 is broad.
Yeah, but you'll need to spend a bit more time with the data to ensure that the peaks are reasonable (you'll have to tweak the parameters if not).

Adding to Devon Ryan's answer, see this image from ENCODE.


text:
Issues performing variant calling with GATK
I am trying to perform variant calling on a BAM file generated through STAR version STAR_2.6.0b for wheat genome using GATK haplotypecaller as follows:
gatk HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar HaplotypeCaller -I sorted.bam -R wheat.fa -O st.vc
02:34:13.450 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
02:34:13.600 INFO  HaplotypeCaller - ------------------------------------------------------------
02:34:13.600 INFO  HaplotypeCaller - The Genome Analysis Toolkit (GATK) v4.0.4.0
02:34:13.600 INFO  HaplotypeCaller - For support and documentation go to https://software.broadinstitute.org/gatk/
02:34:13.601 INFO  HaplotypeCaller - Initializing engine
02:34:40.721 INFO  HaplotypeCaller - Done initializing engine
02:34:48.175 INFO  HaplotypeCallerEngine - Disabling physical phasing, which is supported only for reference-model confidence output
02:34:48.598 INFO  HaplotypeCaller - Shutting down engine
[June 7, 2018 2:34:48 AM PKT] org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 0.59 minutes.
Runtime.totalMemory()=10305929216
java.lang.IllegalArgumentException: samples cannot be empty

When I searched for this error the most proposed solution was to validate BAM file as suggested here. The ValidateSamFile gave following output for the bam file:
gatk ValidateSamFile -I sorted.bam -M SUMMARY
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar ValidateSamFile -I sorted.bam -M SUMMARY
02:53:27.498 INFO  NativeLibraryLoader - Loading libgkl_compression.so   from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 02:53:27 PKT 2018] ValidateSamFile  --INPUT sorted.bam --MODE   SUMMARY  --MAX_OUTPUT 100 --IGNORE_WARNINGS false --VALIDATE_INDEX true --INDEX_VALIDATION_STRINGENCY EXHAUSTIVE --IS_BISULFITE_SEQUENCED false --MAX_OPEN_TEMP_FILES 8000 --SKIP_MATE_VALIDATION false --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 02:54:08     SamFileValidator        Validated Read    10,000,000 records.  Elapsed time: 00:00:36s.  Time for last 10,000,000:   36s.  Last read position: TGACv1_scaffold_080127_1DS:28,522
INFO    2018-06-07 02:54:44     SamFileValidator        Validated Read    20,000,000 records.  Elapsed time: 00:01:12s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_602746_7DL:127,507
INFO    2018-06-07 02:55:20     SamFileValidator        Validated Read    30,000,000 records.  Elapsed time: 00:01:48s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_641946_U:17,161
INFO    2018-06-07 02:55:55     SamFileValidator        Validated Read    40,000,000 records.  Elapsed time: 00:02:23s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_471755_6AL:43,097
INFO    2018-06-07 02:56:31     SamFileValidator        Validated Read    50,000,000 records.  Elapsed time: 00:02:59s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_434075_5DL:29,039
INFO    2018-06-07 02:57:07     SamFileValidator        Validated Read    60,000,000 records.  Elapsed time: 00:03:34s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_062491_1DL:40,466
INFO    2018-06-07 02:57:42     SamFileValidator        Validated Read    70,000,000 records.  Elapsed time: 00:04:10s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_558158_7AL:13,072
INFO    2018-06-07 02:58:18     SamFileValidator        Validated Read    80,000,000 records.  Elapsed time: 00:04:46s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_435811_5DL:29,008
INFO    2018-06-07 02:58:53     SamFileValidator        Validated Read    90,000,000 records.  Elapsed time: 00:05:21s.  Time for last 10,000,000:   35s.  Last read position: TGACv1_scaffold_572030_7AS:9,114

## HISTOGRAM    java.lang.String
Error Type      Count
ERROR:MATE_NOT_FOUND    107771
ERROR:MISSING_READ_GROUP        1
WARNING:MISSING_TAG_NM  96856336
WARNING:RECORD_MISSING_READ_GROUP       96856336

In order to remove the errors shown by ValidateSamFile I used FixMateInformation as suggested here. as shown below(sharing part of file):
gatk FixMateInformation -I sorted.bam -O mate.bam -MC
Using GATK jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar
Running:
java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /home/geslab/bin/gatk-package-4.0.4.0-local.jar FixMateInformation -I sorted.bam -O mate.bam -MC
06:13:13.934 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/home/geslab/bin/gatk-package-4.0.4.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
[Thu Jun 07 06:13:19 PKT 2018] FixMateInformation  --INPUT sorted.bam --OUTPUT mate.bam --ADD_MATE_CIGAR true  --ASSUME_SORTED false --IGNORE_MISSING_MATES true --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 500000 --CREATE_INDEX false --CREATE_MD5_FILE false --GA4GH_CLIENT_SECRETS client_secrets.json --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false
INFO    2018-06-07 06:13:24     FixMateInformation      Sorting input into queryname order.
INFO    2018-06-07 06:23:26     FixMateInformation      Sorting by queryname complete.
INFO    2018-06-07 06:23:26     FixMateInformation      Output will be sorted by coordinate
INFO    2018-06-07 06:23:26     FixMateInformation      Traversing query name sorted records and fixing up mate pair information.
INFO    2018-06-07 06:23:33     FixMateInformation      Processed      1,000,000 records.  Elapsed time: 00:00:06s.  Time for last 1,000,000:    6s.    Last read position: TGACv1_scaffold_361122_4DS:70,952
INFO    2018-06-07 06:23:42     FixMateInformation      Processed     2,000,000 records.  Elapsed time: 00:00:15s.  Time for last 1,000,000:    8s.    Last read position: TGACv1_scaffold_044571_1BL:711
INFO    2018-06-07 06:23:50     FixMateInformation      Processed     3,000,000 records.  Elapsed time: 00:00:23s.  Time for last 1,000,000:    8s.   Last read position: TGACv1_scaffold_433135_5DL:36,647
 .
 .
 .
 .
 [Thu Jun 07 06:40:50 PKT 2018] picard.sam.FixMateInformation done.    Elapsed time: 27.62 minutes.
 Runtime.totalMemory()=8004304896
 Tool returned:
 0

After performing FixMateInformation I again validated the BAM file using ValidateSamFile but the errors are still the same as shown above which shows Mate not found, which means error is not fixed.
How these errors can be fixed? So that I can perform variant calling using BAM file?

The primary error isn't the missing mate, but likely the missing read groups. Use Picard AddOrReplaceReadGroups to ensure your BAM files have read groups and then run those through the haplotype caller.

Should I run AddOrReplaceReadGroups with default parameters as discussed here or there is something that needs to be extracted and provided as parameters for AddOrReplaceReadGroups?

There are no defaults for the read group information. You can just make up values to put there (usually I just use the sample name). Note that you only need the required fields.

You could add read groups:
['STAR', '--twopassMode', 'Basic', '--genomeDir', '/data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/',
         '--runThreadN', '10', '--limitSjdbInsertNsj', '5000000', '--outSAMtype', 'BAM', 'SortedByCoordinate', '--twopass1readsN', '-1',
         '--sjdbOverhang', '100', '--readFilesIn', sra1 + '_1.fastq', sra1 + '_2.fastq',
         '--outSAMattrRGline', 'ID:' + rg1, 'SM:' + rg1, 'PL:ILLUMINA'], shell=False)

In other words, "add read groups".

@user3018491 on which version of STAR you have used the parameter SortedByCoordinate? As discussed here I tried to use this on version 2.6.0b but it failed.

STAR
Usage: STAR  [options]... --genomeDir REFERENCE   --readFilesIn R1.fq R2.fq
Spliced Transcripts Alignment to a Reference (c) Alexander Dobin, 2009-2015
versions
versionSTAR             020201

The problem is that you are missing read groups. You can add them directly to bam files you already have using picard tools. Something like
java -jar ./software/picard.jar AddOrReplaceReadGroups \
  I=input.bam \
  O=output.bam \
  RGLB=lib1 \
  RGPL=illumina \
  RGPU=unit1 \
  RGSM=hello

Hello 罗素幸福之路, it would be nice to explain what the code you have posted means. I tried to add at least one sentence, feel free to [edit] it and improve it further.


text:
Seurat CCA, AlignSubspace, Argument dist.method is only useful with multivariate timeseries
I am trying to combine two datasets to be able to compare the cells in them and using Seurat CCA for it:
https://satijalab.org/seurat/Seurat_AlignmentTutorial.html
When I am running AlignSubspace function I am getting a warning for each single dimension that I selected it to be run on:
seurat_union <- AlignSubspace(object = seurat_union, reduction.type =
    "cca", grouping.var = "protocol", dims.align = 1:13)

So, I am getting 13 warnings:
> warnings()
Warning messages:
1: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
2: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
3: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
4: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
5: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
6: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
7: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
8: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
9: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
10: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
11: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
12: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
13: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries

What is this warning telling me? What should I pay attention to with my datasets?

Could you please clarify if you had these warnings running the tutorial or with your own dataset? If it is the second case could you please provide information about your data if it is a time series or not, and the type of data you add as an input?

With my own dataset. The input is two previously saved Seurat objects. I did analysis of the datasets separately, but then needed to combine them to compare. The input to Seurat objects were raw expression matrices for single cell rna seq data.

And are you sure that seurat_union is the type of the object it is needed. It mind if it wouldn't be you would get an error, but a warning might indicate that perhaps you need to add some other type of data with this object or with other attributes. But it simply might be that the dims.align parameter should be something else

I am following the tutorial step by step except that I have already treated the input Seurat objects: normalized them, did scaling, clustering, tsne, etc. they way it is specified in a different Seurat tutorial.

I am able to align them, and work further, it does not give any other error down the line, but I am afraid that something still went wrong, the alignment was not as good as expected and maybe I will get some result at the end but it may be flawed because I did not take into account this warning - that is what I am afraid of

I don't have enough reputation to comment so I will briefly answer here
This warning is fixed in Seurat 2.3.2. Without upgrading you can ignore this message. It doesn't mean anything for the actual analysis performed so your results won't be flawed.


text:
Using R to import specific records from a PDB file
I'm using R and I'm trying to import specific records from a PDB file (e.g. HETNAM, HETSYN). Unfortunately, it seems the read.pdb function of either biomaRt or Rpdb aren't capable of it.
Any ideas?
Thanks

Hi, could you show what have you tried? As far as I know it is not so difficult to parse the files, but that depends on the format of the file you downloaded.

Hi, nothing special, I just used the read.pdb function to read in the PDB files.

If your file is in PDB format as the one shown here in .pdb format then you can simply extract your desired fields using sed which is a programming language used to parse and transform text. As for this structure i.e 6C4A in .pdb format HETNAM and HETSYN can be extracted using sed as follows:
sed  -n -e '/HETNAM/ p' -e '/HETSYN/ p'  6C4A.pdb

HETNAM     EJA S-[(1Z)-2-CARBOXY-N-HYDROXYETHANIMIDOYL]-L-CYSTEINE
HETNAM     PYR PYRUVIC ACID
HETNAM      MG MAGNESIUM ION
HETNAM     3NP 3-NITROPROPANOIC ACID
HETNAM     PEG DI(HYDROXYETHYL)ETHER
HETNAM     EPE 4-(2-HYDROXYETHYL)-1-PIPERAZINE ETHANESULFONIC ACID
HETSYN     EPE HEPES

Thanks, seems pretty straightforward. How do I use sed? Is it a way to apply it on R ?

It is straight forward. Write this command on Linux command line/terminal with your input file present in current working directory.

Secondly for using sed in R you should check this