text:
What value of reference RMSD after docking simulation from any starting coordinate would be considered outrageous to have?
I need to lightly verify a set of docking simulations from the reference RMSD of the outputs, and as we all know reference RMSD depends on the initial coordinates of the ligand set. What value of reference RMSD after docking from any starting coordinate would be considered outrageous to have?

Based on your terminology I am guessing you are using Autodock Vina. Right?

@MatteoFerla Any docking software actually, but do u have an idea for Autodock Vina? I am trying to find a general RMSD pose range for results in docking simulation

I think I follow your question now, sorry. It utterly depends how much you know of where the ligand actually goes.
I don't like VINA and use Rosetta would use the score command with the reference structure given with the -in:file:native flag, while all the other poses as -s * as normal. But in a Rosetta docking_protocol run how much a ligand wonders off course depends on the flag you set, hence my confusion.

Let's get some definitions out of the way:
RMSD, root mean square deviation, is a metric of distance of molecule A and molecule B. Think of it as average euclidean distance. It us casually called just 'deviation', but has nothing to do with standard deviation (the root of variance, 2nd moment). The units are Ångström.
Reference RMSD is a term that appears in some manuals and is a the RMSD between the calculated pose and the starting pose. When you set up your docking you as a user define the starting position and how much it can deviate from that position. If you have a crystal structure with a ligand and you are docking the same or similar compound, then yes, it makes total sense to worry about this metric. In which case, RMSD greater than 2 are worrisome in a 2 Å resolution structure. You especially see RMSD as a metric of accuracy when PDB structures are redocked, wherein the same ligand is docked and reference is the file from the PDB (or its PDB Redo derivative), so RMSD greater than 2 Å is bad as it would not match the density, while under 1 is really good —note that in tests about 50% and 75% of re-docked ligands are docked less than 2 with AutoDock (free) and ICM-Dock (pay) respectively. But if you are doing it for an experiment and not for a test, you would probably have set the docking to be highly constrained in the first place (the grid size in AutoDock or position randomisation factor in most others, eg. Rosetta Ligand dock). If you actually care about reference RMSD the thing to manually look at is the top 3-4 poses in the electron density of your native structure. If you don't really know how your ligand docks please ignore this value as it is just a product of your starting pose. It is not a measure of how "happy" the complex is.
ΔGibbs is a negative value in kcal/mol that is basically a metric of how "happy" the atoms are. It is actually an energy potential, but confer Wikipedia pages for the proper definition. The more negative it is the happier the atoms are. Force-field calculations solve this value and the difference between undocked and docked (ΔΔG) is a good metric of how well your ligand docks. Please use this.
Affinity ΔΔG is not a k_off or a k_on. That is a result of the ligand binding to nearby positions and how well it can get into the binding site (e.g. in P450 enzymes the ligand orientation is dictated by how it squeezes into the entrance way). And how insoluble the ligand is, which is a major problem*. For a better metric, there is dynamic undocking where a constrained MD simulation is done with the ligand getting pulled out of the pocket. But this is not trivial.
∗ As any comp-chemist can tell you it is mighty hard to calculate solubility. Some force-field calculations do not even take account of π-π or S-π interactions. Drug-like compounds are rigid aromatic and rigid and highly hydrophobic. Many programs use a hack called implicit waters, where the water is not simulated as individual molecules but as a flowing thing. This means that ΔΔG benefits of constrained waters are ignored and alternative protontation states of water too. This is fine for protein (mostly), but has major errors for drugs. But in terms of thing to look out for in a docking algorithm, explicit waters is not the first as conformer sampling and side chain perturbation are more important.

Thanks for your response

2A would be huge for docking, but agreed depends on the resolution.
Its 10A if you are looking at protein-protein bonds in a beta pleated sheet and your molecules need to be a lot smaller than that.


text:
How to edit the headers of multiple fasta files from multiple folders? (recursively)
This is an expansion of my previous question, How to edit the headers of multiple fasta files from multiple folders?
My directories are organized as follows: one main directory, in which I have multiple project directories. Inside each project directory, there are multiple metagenome directories. Inside each metagenome directory, there is always a single .fna (fasta) file.
Names: project directories start with 'mgp', metagenome directories start with 'mgm' and end with '.3'. They end in numbers, but the numbers are not related. (eg. project directory 'mgp83581' contains metagenome directory 'mgm4729322.3').
What I want is a loop that would go through all these fasta files, and add to the headers codes corresponding to project and metagenome directory. For example, if the fasta file in project directory mgp83581, metagenome subdirectory mgm4729322.3 had the following headers...
>seq1
>seq2
...then I would want to change it to:
>83581_322_seq1
>83581_322_seq2
Basically, I want to add the number code associated with the project, and the last three digits of the 7-digit code associated with the metagenome (apart from the '.3' part of the file name).
Below is the answer that I got for my previous question, where I didn't have project directories, and where the 3-digit code for the metagenomes stayed in a manageable range. (thanks terdon!)
for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done
Using ubuntu. I would appreciate your help!

You can edit your for loop so as to include the paths to your files of interest.
The command below will give you the relative paths (relative to where you run the command, I think it is better to run this within your main directory) of your files of interest:
ls mpg*/mgm*.3/*.fna

You can use the output of this command (your files of interest) to iterate over via your for loop. Your code will look like:
for i in {...}; do
    sed -i.bak "s/^>/>${i}_/" i
done

Alot slicker than using Perl

Thank you! I've tried to do this. (though tried 'while' instead)
ls mgp*/mgm*.3/*.fna | while read line ;  do  header=$(sed -r "s/^mgp([0-9]*)\/mgm[0-9]{4}([0-9]{1,})\.3\/(mgm.*\.fna)/\1_\2/")   sed -i.bak "s/^>/>${header}_/" $line done
It seems to only edit the first of the files I'm interested in (but it does it correctly as far as I can tell). Any idea of what could've gone wrong?

I am not sure if you really need that complicated command within do, fasta headers always start with >, which seems to be parsed in your sed. My suggestion would be: ls mgp*/mgm*.3/*.fna | while read line; do sed -i.bak "s/^>/>${i}_/" $line; done.

I've tried. For some reason it doesn't work for me. All it does is add an underscore to the fasta headers in the edited files.

It is because I forgot to change the  ${i} part in sed while editing the code regarding your comment above with while read. You can do something like ${line:Ind:Len}, Ind being the zero based index and Len is the length of the substring you would like to extract from the full file path.

source: Stack Exchange

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

sort: votes


text:
Merging two files together
Is it possible to merge two files together using awk? I have two files with matched  Chrom and pos  column
file1 has the columns
CHROM   POS ID  REF ALT QUAL    FILTER  INFO

file2 has the columns
CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

I want the final merged file to have these columns
CHROM   POS ID  REF ALT QUAL    FILTER  INFO FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

Can you show us a few more lines from each file? Do the files have other headers or is it only the fields? If they're real vcf files, they will have headers which complicate things. If they're headerless lists of variants, there are better tools than awk.

Matched CHROM and POS alone might create problems if there are multi-allelic variants (or even if ALT is different between the two files). You might want to match all three of those columns to be sure that the variants you're merging are identical.

If you just need to merge column together, you can use paste command:
paste file1 <(cut -d\  -f6- file2)

In this case the whole file1 is taken and only part of file2 starting from column 6.
Assuming your data are sorted, another alternative is to use join:
join -o '1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 2.6 2.7 2.8 2.9 2.10 2.11 2.12' file1 file2

The -o option allows to select column from either file1 or file2.

I'd recommend using R and either merge or one of the joins. Assuming your data is not too large, data.frame + read.delim should suffice to hold the two datasets. If not, try data.table + fread..

I'd agree with this.  R merge will work even if things are out of order, or some things are missing from one file or another.

I did something like  below,
awk 'FNR==NR{a[$1$2]=$8 FS $12;next}{ print $0, a[$1 $2]}' file1.vcf file2.vcf > Final.vcf

But then the output looks like
CHROM   POS ID  REF ALT QUAL    FILTER  FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2 INFO

is it possible to arrange the output to look like
CHROM   POS ID  REF ALT QUAL    FILTER  INFO FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

It's not awk: but this is a fairly simple python/pandas script (assuming tab-delimited).
file1
CHROM   POS ID  REF ALT QUAL    FILTER  INFO

file2
CHROM POS ID REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2

#!/usr/bin/python
import pandas as pd

file1 = pd.read_csv('file1.txt', sep='\t')
file2 = pd.read_csv('file2.txt', sep='\t')

# you can compare file1.shape vs. file2.shape to ensure they match
# you can also use fileN.sort_values(by=['POS', 'CHROM'], axis=1) to ensure sort matches

file3 = pd.concat([file1, file2[['FATHER', 'MOTHER', 'DAUGHTER1', 'DAUGHTER2', 'DAUGHTER3', 'SON1', 'SON2']]], axis=1)

file3.to_csv('merged_file.txt', sep='\t', index=False)

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


text:
How to troubleshoot access to tabix file via http or https?
I have some BED files that I have compressed with bgzip and indexed with tabix.
I am using version 1.9-392-gcfa32f7-dirty on the server side.
$ /usr/local/bin/bgzip --version
bgzip (htslib) 1.9-392-gcfa32f7-dirty
Copyright (C) 2019 Genome Research Ltd.

The files are served from a subdirectory made available via nginx on port 80 (http) and 443 (https).
The problem is that, for some bgzip-compressed files, I get a correct response to a tabix query if I run it locally, e.g.:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KL.gz" chr1:50669400-50669800
chr1    50669200    50669400    id:253347,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]
chr1    50669400    50669600    id:253348,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]
chr1    50669600    50669800    id:253349,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]

However, when I run this command from another host (such as my laptop), I get no result:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KL.gz" chr1:50669400-50669800
$

Running echo $? gives a result of 0 (no error). There is nothing in nginx error logs to suggest an issue with http access.
Other bgzip-compressed files work, which further suggests a problem with either the compressed file or its tbi index file, and not with the nginx configuration.
From my laptop, for instance, the following file can be accessed:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KLs.gz" chr1:50669400-50669800
chr1    50669200    50669400    id:253347,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]
chr1    50669400    50669600    id:253348,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]
chr1    50669600    50669800    id:253349,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]

I've investigated the two test files and they are both sorted lexicographically by chromosome name and numerically by start and stop coordinates. I do not see any odd line endings. These two test files are both virtually identical, except for the data in the fourth column.
Is there anything I can do to further troubleshoot this and figure out why I cannot get a result from other hosts other than the server itself, and only for certain files? (Is there a debug version of htslib applications that would include verbose logging, say?)

Update
In response to Jonathan Moore's comment, I wget the bgzip-ed gz archive and its associated tbi index, and then try to query the file locally (not via http, but directly via the filesystem).
I get the following error message:
$ tabix hg19.15.all.KL.gz chr1:50669400-50669800
[W::hts_idx_load2] The index file is older than the data file: hg19.15.all.KL.gz.tbi

The tbi index file was made immediately after creating the gz file, so I don't understand this error.

Can you wget these particular files from the laptop?  If so, can you then use tabix locally on the laptop to decode the file you just got?

When I do so, I get a strange error. I have updated my question with more details.

That looks to be a fairly common warning message that people get in a few different strange circumstances.  The only thing that comes to mind is whether there is a weird version difference between the tabix format installed on the laptop and that on the server that is coming into play.  Can you get the tabix's to tell you their version numbers?

My other thought is, if you generate a tabix index on the laptop, is it different from the one that you wget from the server?


text:
match two files and check if position of chrom is within exon region
I have two files, file1 and file2. I want to match file1 chrom with file2 chrom, and if there is a match then I want to check file2 if the position of the match chrom1 is within the exon region that is in file1. i want to do this with the whole file and print out rows that match and the position is within the exonic regions.
For instance, file1 has a chrom 1 and POS 736689 which match chrom 1 in file2, now I take POS 736689 from file1 and check if it is within the exon region 29553,30563,30975,   30039,30667,31097,... in file2.
If it is within the region it will print out the file1 row that matches the pattern.
The exon regions can be read as:
29553 to 30039
30563 to 30667
30975 to 31097

file1 look like this:
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO
1   736689  T   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=6
1   758663  CGCC    TTAG    0/0 0/1 0|0 0|1 0|1 0|0 0|1 AC=2;AN=2184
1   762005  G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=13
1   771950  C   T   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=110;AN=2184
1   772172  C   T   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=8
1   790590  C   T   0/1 0/1 0|0 0/1 0/1 0|0 1|1 AN=2184;AC=84

file2 look like this
chrom   Start   End
chr1    29553,30563,30975,  30039,30667,31097,
chr1    53048,54829,    53067,54936,
chr1    69090,  70008,
chr1    160445,161313,  160690,161525,
chr1    324755,327551,  326514,328453,
chr1    367639, 368634,
chr1    459655,461153,461750,   459992,461288,461954,
chr1    523008,529614,  523081,530148,

Please edit your question so that your data is properly formatted (like you did withh file headers) and provide some easily parsable example data. I have posted an answer but not all people would (understandably) spend time cleaning data before suggesting a solution.

Despite an awk solution is favored by the OP, here is an R solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward awk solution, at least for me).
res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84

File 1 looks like a vcf file? Maybe you can read it in using VariantAnnotation package from Bioconductor, and do overlap with the second file. To make it easier, one can convert the second file to a GRanges object in R.

OP seems to add awk tags to all their posts even though multiple people have mentioned that R would be a better way to merge/lookup files.

It's a good answer by @haci nevertheless

BEDTools already has a function to take a vcf and bedfile as input and look for intersections.  Most people would get the exon coordinates into bed format, and use that.

Make a sorted BED file called exons.bed and run BEDOPS bedmap to map variants to exons:
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed <(vcf2bed < variants.vcf) > answer.bed

To make a sorted BED file from your exons, I'd suggest using sort-bed.
I would also suggest using awk to make sure that the chromosome naming scheme is consistent between exons and variants, e.g.:
$ vcf2bed < variants.vcf | awk '{ print 'chr$0' }' > variants.bed

Then use that result, instead, for mapping:
$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed variants.bed > answer.bed

If chromosome names are not consistent, mapping or other set operations with BEDOPS or other toolkits will not work.


text:
Selecting genes with more contribution from PCA
I have RNA-seq data in response to treatment vs non response; By machine learning I selected three principle components likely can predict the response based on the gene expression. Now I have something like this
> head(pca$rotation[1:10,1:6])
               PC1          PC2          PC3          PC4           PC5          PC6
ACTB   -0.01858031 -0.005810636  0.014939832 -0.009594126 -0.0118313644 -0.011534401
ATP5F1 -0.01822611  0.014851900 -0.016401920 -0.015616379  0.0254850111 -0.010225516
DDX5   -0.02477827 -0.020573978  0.002199966 -0.009296852  0.0006222345  0.014506614
EEF1G  -0.02803284  0.001126758 -0.010570521 -0.017035205  0.0007583815  0.011225639
GAPDH  -0.02691309  0.005503228  0.006446502 -0.012686738 -0.0289737909  0.007036631
NCL    -0.03019765 -0.003011778 -0.015534611 -0.002019556 -0.0126401245  0.013977142
>
> dim(pca$rotation)
[1] 3962   56
>

Machine learning says PC 6, PC 21 and PC 39 are better; how I know from 3962 genes in each of these PCs which genes is more contribute to response to treatment?

Have you looked at the loadings in those PCs?

Hi @Fereshteh you should know this, you have enough experience in ML. PCA cannot solve this but you can use the classification for supervised learning which will tell you

Thank you, yes; But I can not conclude which genes I should select from preferred PCs

Like @DevonRyan mentioned, you can go by the loadings, which will be sort( abs(pca$rotation[,"PC6"]),decreasing=T) for example in your case. But it kind of defeats the purpose of PCA right? You use PCA to get linear combinations of the genes..

Thank you @StupidWolf; So you think genes with higher loading scores in each PCs don't necessarily can predict the response to treatment as here I am supposing a linear relationships between genes?

Yes the genes with the higher loadings in the better PCs should contribute more towards response / non-response. My point is, it should be a combination of the genes. For example, if PC6 has strong +ve loadings for gene 1,2,3 and strong -ve loadings for gene 4,5,6, I would interpret it as, a strong expression of genes 1-3 and low expression of genes 4-6 contribute towards a response phenotype

What exactly is the goal of your analysis? If it is to predict and then point out a few potential causative genes, then the PCA approach you have is ok.
I would be cautious about using up to PC39 (which explains very little of the variance in gene expression).

Thank you @StupidWolf; Yes as you mentioned the ultimate goal is point out a few potential causative genes that predict is a patient would respond to treatment or not. I have 56 PCs (I have 56 samples in total). I have three PCs that gives a final models with the highest accuracy then I need to know more important genes underlying these PCs. Here by using prcomp I have some PCs then I used these PCs as independent variable to predict response to treatment finally three PCs showed can predict the responder patients better.


text:
How to define an outgroup to build a robust amino acid tree
I building a robust amino acid phylogeny with adequate robustness support (see previous post here).
This is a brief summary of what I have done:

I performed a BLASTp analysis using a protein sequence from species A as query against the full NCBI protein database,
downloaded all the hits in fasta format and I added the sequences I identified with the query sequence in my target genome, which is closely related to species A.
I then aligned the sequences using the Compute - Alignment - MUSCLE option in Seqotron and I saved my alignment file in Phylip format.
I tried to run this command line in RAxML to build a robust maximum likelihood tree:

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s ~/Desktop/alignment_file.phy
However, the program showed the following error:
Error, the outgroup name "outgroup-p" you specified can not be found in the alignment, exiting ....
Question how can I build an outgroup for my phylogenetic tree, and how do I include it in my alignment?
I appreciate your help and suggestions!

Can you show us the input file please? Personally I would avoid using a '-' in the filename

Sure. The input file is named simply "alignment.phy" to avoid any problems.

I think you have a bug
raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /Users)/username/Desktop/alignment_file.phy

Try the above, where username is your username. I'll check my codes later. The most important bit is the space between the -p and word outgroup. You can omit the -p 10000 it is only needed  for very geeky maximum likelihood and sets the random number stream.

Thank you @Michael G. I tried using this line:  raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /Users/username/Desktop/alignment_file.phy but I get the following error: The file /Users/username/Desktop/alignment_file.phy you want to open for reading does not exist, exiting ...  I also tried using raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /username/Desktop/alignment_file.phy but I got the same error message. As mentioned above, the file name is simply alignment.phy.

I've set up a chat group here, https://chat.stackexchange.com/rooms/99546/phylogeny. Please drop by and let me know when you are around. Tomorrow afternoon your time zone would be good. Please make sure you have access to your Linux/Darwin terminal and your data sets

Thanks a lot @Michael G. I will be there on Tuesday October 8 from 17:00 (CDT) onward. Does that work for you? I can connect earlier if it suits you better.

Yes earlier, ie before 17:00 please (6 hour time difference). Must leave by 17:00

I see. How about 9:30 am CDT (UTC -5), which would be 15:30 BST (UTC +1)? Is this better?

Ok, see you then

To clear up what may be a possible misunderstanding about RAxML: RAxML (and most/all other maximum likelihood phylogenetic inference programs using reversible models of sequence evolution) explicitly infer unrooted trees. This is because, for a reversible model of sequence evolution, changing the position of the root does not affect the likelihood of the tree. Therefore, as you will see in the RAxML manual, specifying an outgroup is explicitly a display option - it doesn't affect the inference in any way. This means that it's usually easier to just infer the tree in RAxML without specifying any outgroup, and then use any tree viewing software (e.g. FigTree, phylotree.js, or equivalent) to reroot the tree as you need to.


text:
vcf tools Error: Require Genotypes in VCF file in order to output Fst statistics
I want to calculate Fst by vcf tools and GATK.
I did this steps:

for creating gvcf:

gatk HaplotypeCaller -R ref.scf.fasta -I input.bam -ERC GVCF -O out.g.vcf

for combine:

gatk CombineGVCFs -R ref.scf.fasta --variant 1.g.vcf --variant 2.g.vcf --variant 3.g.vcf --variant 4.g.vcf --variant 5.g.vcf --variant 6.g.vcf --variant 7.g.vcf --variant 8.g.vcf -O outt.vcf

for genotype:

gatk GenotypeGVCFs -R ref.scf.fasta --variant outt -O outt.vcflist

for Fst:

vcftools --vcf outt.vcflist --weir-fst-pop 1.txt --weir-fst-pop 2.txt --out 1_vs_2
There was no error till step 3 and successfully done. But, I got error in step 4:
Error: Require Genotypes in VCF file in order to output Fst statistics.

I used the last version (GATK/4.1.3.0-Java-1.8) and the last version of vcf tools (VCFtools/0.1.16), but I still have the same error!
Would you please advise me?
Best,
Razi


text:
too long header for fasta file
I have a fasta file, like this:
>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039::m.3039 TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3::g.3039  ORF type:complete len:100 (-) TRINITY_DN100000_c1_g1_i3:1027-1326(-)
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365::m.84365 TRINITY_DN100000_c1_g6::TRINITY_DN100000_c1_g6_i2::g.84365  ORF type:complete len:112 (-) TRINITY_DN100000_c1_g6_i2:379-714(-)
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

if I want to only keep the second header, like this (the amino acid sequence keep unchanged). which command could i USE?
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

I know a command
> sed 's/::.*//' input > output

can only keep the first header >TRINITY_DN100000_c1_g1 if I want to keep the second header with the isoform information >TRINITY_DN100000_c1_g1_i3, how should I correct this command?

You can also do this via eFetch or Biopython directly without regex-ing

Here's one way:
$ sed -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

The idea is to match a >, then as many non: as possible ([^:]*), and then capture the string between the first :: and the second ::. Finally, you replace the entire line with the captured string.
Once you're satisfied this does what you need, you can edit the original file with -i:
sed -i -E 's/>[^:]*::([^:]*)::.*/>\1/' file.fa

Alternatively, set awk's field separator to :: and print the second field only on lines that start with >:
$ awk -F'::' '/^>/ ? $0=">"$2 : 1' file.fa
>TRINITY_DN100000_c1_g1_i3
MVWIKFRGLHRVLTSTPLVKSGKTPSQTWAFLDISVELIVFLFLNVHKSPMPHFKIYSEA
FSEEWSLLWLQYSRHLIQKPKPWQIKIELLHLCCCNRLC*
>TRINITY_DN100000_c1_g6_i2
MEMMQEIIPFAREMLSARPSKGTMKVYLVGGTFAVLGIVSGMVEAACSLFPEQEESTLTK
LMEDCLTVTAQNQEPQTFIPEDDEQDAEMEAKAKDLPMFRQRRMSFRAHAS*

A 2 step workaround is:

Get rid of everything after the piece of info you want, that is anything starting with :: and goes on with g:
sed "s/::g.*//" fasta > fasta2

>TRINITY_DN100000_c1_g1::TRINITY_DN100000_c1_g1_i3

Get rid of the "first header" as you call it, after the first step, it is rather straightforward:
sed "s/>TRINITY_.*::/>/" fasta2 > fasta3

>TRINITY_DN100000_c1_g1_i3
EDIT:
@terdon showed that multiple sed commands can be used when separated by ; in GNU sed (default on Linux):
sed "s/::g.*//; s/>TRINITY_.*::/>/" fasta > modified_fasta

it worked!! thank you!

There is no need for multiple steps: sed 's/::g.*//; s/>TRINITY_.*::/>/'. Or, if you don't have GNU sed: sed -e 's/::g.*//' -e 's/>TRINITY_.*::/>/'. You could even pipe one sed into another, that would at least avoid the temp files: sed 's/::g.*//' |   sed 's/>TRINITY_.*::/>/'.


text:
Is it okay to use deeptools bamCompare (SES normalization) for comparisons across different ATAC-Seq datasets?
We are trying to use deeptools for analysis of ATAC Seq datasets. We have datasets with different sequencing depths and are wondering if bamCompare's SES based normalization is appropriate for comparing two ATACSeq datasets as opposed to ChipSeq datasets. If such normalization is okay to perform, does it make sense to compare difference in values across comparisons? Let's say I do bamCompare on A and B using SES, and then I do bamCompare (using SES) on A and C. Will it then make sense to compare the two comparisons (A-B & A-C)?
If such normalization isn't okay, what are the best tools to normalize ATAC-Seq datasets with different sequencing depths so that one can perform comparisons like above?
Thanks for any suggestions or comments.

I would be very hesitant to use SES normalization in a case like this. We recommend that people run plotFingerprint first and see if there is good separation between the samples and only then using SES. I would not expect there to be drastic separation between the samples for ATAC-seq data, so I don't think it will work well. Our standard ATAC-seq normalization is RPGC (aka, "1x normalization"). That should work well for your use-case as well.
As an aside, if you really want to make statistically-backed comparisons between groups I can recommend the CSAW package.

Thanks Devon! multiBamSummary tool produces a scalingFactor table for respective BAMs. Does it make sense to first scale with these scaling factors, create bigwig using bamConverge and then use these bigwigs for comparison using bigwigcompare? I am wondering if scaled files go to bigwigcompare, will it make for a better comparison across comparisons.

Yes, that's another good way of approaching things. bigwigCompare does assume that the input bigWigs are scaled appropriately and the procedure you suggest will create a nicely robust comparison.

Thanks! You mentioned plotFingerprint for seeing good separation. I did the plotCorrelation heatmap. Is there some approximate threshold for "good separation"? Let's say I have a correlation between 0.5 and 0.7 across different samples. Will this make them ok for SES?

There's no substitute for plotFingerprint here.

Devon,
Can we use scalingFactor table (from multiBamSummary) to create new scaled BAMs instead of bigwig/bedgraph? I want to perform peak calling only on normalized/scaled BAMs so I can properly compare peaks/footprints across different samples. I know I can downsample using samtools but some scaling factors are greater than 1 so not sure how to "upscale" it.
Thanks again!

So  I found a post where you suggest not to upsample:
https://www.biostars.org/p/365825/
But lets say after running multiBamSummary, i get a table that looks like this:
A     1.1215
B     1.4657
C      1.3783
D      1.0529
E      0.6131
F    0.5919
G       0.9988
H   1.9739
I   1.0636

In order to make these samples comparable, do you have any suggestions on how do I go about downsampling everything without upsampling?

At least some peak callers can accept bedGraph files, which you can normalize. The issue with normalizing BAM files is that you're quite limited in terms of fidelity (e.g., you can't remove a tenth of a read from a region).


text:
How to load output from the kegg API in biopython into a pandas dataframe?
Biopython provides a (bit unintuitive) API to access to the Kyoto Encyclopedia of Genes and Genomes (KEGG). I am trying to make use of it, but the output is quite unhandy as a string. What is the best way to parse the data into a pandas.DataFrame()?
from Bio.KEGG import REST as Kegg
import io

pd.read_table(kegg.kegg_get('K02545'))

It might be helpful if you provide an example of the output so it is easier understand why this string is not useful, and how it relates to your coded solution.

solution from @sören:
I came up with this solution for my problem:
from Bio.KEGG import REST as kegg

def _get_kegg(kegg_id):
    kegg_output = kegg.kegg_get(kegg_id).read()
    results = {}
    for line in kegg_output.split('\n'):
        splits = line.split()
        if not line.startswith(' '):
            if len(splits) > 0:
                key = splits[0]
                value = ' '.join(splits[1:])
                results[key] = value
        else:
            results[key] += ' '.join(splits)
    return pd.DataFrame(results, index=[kegg_id])

_get_kegg_v = np.vectorize(_get_kegg)

def get_kegg_info(kegg_ids):
    if isinstance(kegg_ids, str):
        kegg_ids = [kegg_ids]
    return pd.concat(_get_kegg_v(kegg_ids), sort=False)


text:
awk find pattern
I am trying to find this pattern on a large file but the output doesn't look good. The output doesn't seem to reflect the pattern.
Father = 1|0 or 1/0

Mother = 110 or 1/0

Daughter 1 = 1|0 or 1/0

Daughter 2 = 1|1 or 1/1

Daughter 3 = 0|0 or 0/0

Son 1 = 1|0 or 1/0

Son 2 = 1|1 or 1/1

The content of the file1 looks like this
CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3

SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   3833321 T   C   1/0 1/0 1|1 1/0 1/0 1|1 0|0 AN=2184;AC=801  0.366758242
1   9009352 A   G   0/1 0/1 1|1 0/1 0/1 1|1 0|0 AC=90;AN=2184   0.041208791
1   9328077 G   A   1/1 0/0 1|0 1|0 1|0 1|0 1|0 AC=1;AN=2184    0.000457875
1   10684412    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=2    0.000915751
1   11782120    G   A   1/0 0/0 0|0 0|0 1|0 1|0 0|0 AC=3;AN=2184    0.001373626
1   12034879    A   T   0/1 1/1 0|1 0|1 0|1 0|1 1|1 AC=3;AN=2184    0.001373626
1   16064465    C   A   1/0 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=4    0.001831502
1   16111170    G   A   1/0 1/0 1/0 0|0 0|0 1/0 1/0 AC=1;AN=2184    0.000457875
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   16258118    C   T   0/1 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=1    0.000457875
1   17966719    A   G   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=0;AN=2184    0
1   19751079    CG  C   0/1 0/1 1|1 0|0 0|0 0/1 0/1 AN=2184;AC=1    0.000457875
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
1   21060356    G   A   1/0 1/0 1/0 1/0 1/0 1|1 0|0 AN=2184;AC=541  0.247710623

My awk script look like this
awk '($8~/1\/1/ || $8~/1\\|1/) && ($9~/0\/0/ ||$9~/0\\|0/) && ($7~/^1\\|0/||$7~/1\/0/) && ($10~/1\\|0/||$7~/1\/0/) && ($10~/0\/0/ ||  $10~/0\\|0/) && ($5~/1\/0/|| $5~/1\\|0/) && ($11~/0\/1/||$11~/0\\|1/)&& ($6~/1\\|0/||$6~/1\/0/) && ($12 ~/1\/1/ ||12~/1\\|1/){print}' file1 >file2

What's wrong with the output? "Doesn't look good" isn't very informative. Please [edit] your question and highlight what is going wrong.

The output doesn't reflect the pattern.

Well yes, I imagine so. But since this is a really complex pattern, it would be very helpful if you could show the output you get and explain why it is wrong. I wrote you the script you're using, and I still don't see what the issue is.

Example output of what you got and what you expected would be helpful, without these one can only guess.
First of all your input does not have a row satisfying the criteria you have given, so I added two dummy rows, one with / and one with | used as separator.
Here is the code I used:
awk '($5~/1\|0/ || $5~/1\/0/) && ($6~/1\|0/ || $6~/1\/0/) && ($7~/1\|0/ || $7~/1\/0/) && ($8~/1\|1/ || $8~/1\/1/) && ($9~/0\|0/ || $9~/0\/0/) && ($10~/1\|0/ || $10~/1\/0/) && ($11~/1\|1/ || $11~/1\/1/) {print}' test.txt

And its output (input is yours + 2 dummy rows):
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy

The code is basically yours, except for:

I escaped | with one \ instead of two
I changed this part: ($10~/1\\|0/||$7~/1\/0/). I think $7 should have been $10
I removed a ^ from one of your conditions, I believe they all should start with ^ though, per @terdon's answer on a related question

I believe the code above does what you want, but when I run code above with two backslashes to escape I get the following output and could not figure out why!
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755
1   16200229    C   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AC=9;AN=2184    0.004120879
1   19955045    G   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=4    0.001831502
Dummy  dummy D D 1/0 1/0 1/0 1/1 0/0 1/0 1/1 dummy dummy
Dummy  dummy D D 1|0 1|0 1|0 1|1 0|0 1|0 1|1 dummy dummy


text:
where to find pathogenicity metadata for E.coli genomes
I have a list of about 900 E.coli genome ids (Genbank plus NCTC ids), e.g.:
NCTC8196_1
NCTC11105_1
CP027464.1
CP027459.1
CP027363.1
CP027388.1
CP027355.1
CP027675.1
CP027445.1
CP027325.1
CP027219.1

and I would like to find pathogenicity metadata for them so that I can classify the molecular data I found with the metadata.
Is there a way to get a table of metadata from these genomes that I can use to associate with my molecular data?
NOTES:
After some searching online, I found Enterobase (http://enterobase.warwick.ac.uk), but it seems the site is mostly dedicated to MLST typing, thus they use different naming conventions than the Genbank id. Maybe there is a way to link the info in Enterobase with Genbank ids somehow?
It seems like Enterobase contains project ids, e.g. PRJNA218110. These project ids can also be found on the other side in the GenBank entries, for which we could associate the GenBank chromosomal ids.
Searching for PRJNA218110 in Enterobase gives a very long list of entries that seem to be from Illumina sequencing data, but in this list, the 'Path/Nonpath' column is always empty.
Enterobase seems to contain lots of entries with useful information, such as: the strain is pathogenic, and produced diarrhea in humans. So a table with entries like the one below would be very useful:
Genbank_id: CPXXXXXXX.1
Path/Nonpath: Pathogenic
Disease: Diarrhea

Genbank_id: CPYYYYYYY.1
Path/Nonpath: Nonpathogenic
Disease: NA

What kind of data, exactly? CP027459.1 is an entire genome. Are you trying to see if this strain is a pathogen? If so, for what species? Are you instead trying to see if mutations in this genome can be pathogenic for the bacterium? Please [edit] and clarify what you're after.

I added a comment now referring to what I can find in Enterobase, thanks @terdon .

Thanks! So you are looking for information on which of these bacteria are pathogens in humans? Humans only or do you also want pathogenicity information on any species?

Humans mainly, but any would be great

you can use BacDive. This is the largest database for bacterial metadata.:
https://bacdive.dsmz.de/
e.g. you can search for the sequence accession "NCTC8196" and you would get the metadata of the according strain:
https://bacdive.dsmz.de/search?search=NCTC8196
What kind of data is available, is highly dependent on the strain. Usually type strains are more rich on data that others.
For further questions don't hesitate to contact us:
contact@bacdive.de


text:
Convert PAF format to SAM/BAM format
I have a bunch of PAF files resulting from the alignments of fastq files on a reference genome with minimap2. I would like to convert them into SAM/BAM format so I can use samtools mpileup to recover the read coverage at filtered positions
So far I've looked at paftools.js and found sam2paf but no paf2sam
I am open to suggestions but I'd like to avoid mapping the reads again.
Edit
I have opened an issue on minimap2's github
Answer from Heng Li :

You need to rerun alignment. No tools exist to convert paf to sam.
why there is no paf2sam in paftools.js

It is technically possible, but someone needs to write the code.

I'm not sure the question is answerable in its current format, there doesn't appear to be enough background information and this often results in substancial misunderstanding between OP and the reply. Lets see what the other think.

Thanks for the feedback. So, as far I as know mpileup can take an alignment file and a list of genomic positions (eg BED) and return statistics for each positions, such as the read coverage. mpileup doesn't work with PAF file so I am looking to convert PAF format to SAM/BAM format. If there are alternatives to mpileup, I am open to the suggestions. Is this better ?

It's really weird that they provide sam2paf but not paf2sam. Could it be that the conversion in that direction isn't possible? Perhaps paf is missing information that would be essential for sam?

@terdon I am not a specialist but it's stated that CIGAR can take different format depending on the input/output options. It might also be not written in the PAF in some case. Is that CIGAR essential to sam ?

I would take Heng Li's advice.  minimap2 is fast. Just add the -a option when you re-run, which outputs the alignment in sam instead of paf

The paf format does not have the CIGAR string and base quality information that are typically present in a BAM/SAM file. While it is in principle possible to write a paf2sam converter, it's going to be very difficult to regenerate the CIGAR string with all the indel information (and maybe it's probably impossible to do so) that you'll need in a BAM/SAM file. If your motivation is to get coverage information, it might make more sense to convert the paf file into a bed file, before using a tool like bedtools to convert it into a coverage track.

note that PAF files can have a "cg:" tag that is a cigar string. it is an optional output from minimap2 for example (there is also cs: tag that is a different encoding of the sequence differences, but cg may be more common)


text:
Convert or extract Seurat object as the input of FateID
I am currently using Seurat to do the scRNA-seq clustering analysis. After this, I am planning to use FateID to do the downstream analysis, could someone teach me how to convert or extract a Seurat object as the input for FateID?
Thanks in advance.

Seurat objects have various slots where different types and forms of data are stored. For FateID you would need an expression matrix (with cells as columns and genes as rows) and this can be extracted via:
seurat_object@assays$RNA@data

Here RNA is the assay (you can have multiple assays per sample, RNA, protein, ...) and data is normalized data. You can extract counts or scaled data instead.
If you are using Seurat v2, drop the @assays$RNA part.

Dear haci, thx, I will try.

Thank you @haci, I have figured out it by myself.


text:
R WGCNA error code
I am using WGCNA to identify Consensus modules. I get a warning message when completing the following step:

Call the network topology analysis function for each set in turn

Warning:
Warning message: executing %dopar% sequentially: no parallel backend registered

Does anyone no what this error code means as I cannot plot the scale free analysis, with soft thresholding.
Thank you
powers = c(seq(4,10, by=1), seq(12,20, by=2))
powerTables = vector(mode= "list", length = nSets)
for( set in 1:nSets)
powerTables[[set]] = list(data=pickSoftThreshold(multiExpr[[set]]$data, powerVector= powers,
verbose =2)[[2]])

Please add your code giving the error and especially the packages you have loaded.

I think you forgot to load (or install) the doParallel package as the error warning message is complaining about the lack of a parallel backend. library("doParallel") should do the trick.

Thank you, I have tried this.. but still getting same error message. Packages I installed now are; edgeR, limma, doParallel, WGCNA

What you see is a warning, not an error. Your calculation will run fine, just slower. Unless you see other errors, you should be able to complete all steps of the analysis.


text:
How to retrieve values from one table that corresponds to the values in another table in R
I have following table.
structure(list(n = c(29L, 11L, 10L, 10L), domain = c(1L, 5L,
6L, 32L)), row.names = c(NA, -4L), class = c("tbl_df", "tbl",
"data.frame"))

I need to retrieve rows from another table that corresponds to the values of the 1st column of the above table.
I did use following code by manually entering the column values of the table.
target <- c("1", "5", "6", "32")
selected_domains = filter(gene_domains, domain %in% target)

gene_domains is the second table.
I have other tables that have many rows and entering manually is impossible. So I will be grateful if someone can help me with this issue.
Thank you.
If you need, this is the second table (gene_domains).
structure(list(gene = c("ENSG00000230594 ENSG00000171155 ENSG00000224089 ENSG00000230347 ENSG00000236446 ENSG00000186471 ENSG00000101892 ENSG00000182890 ENSG00000232119 ENSG00000131721 ENSG00000101882 ENSG00000101883 ENSG00000242362 ENSG00000226685 ENSG00000125352 ENSG00000236126 ENSG00000237957 ENSG00000005893 ENSG00000125355 ENSG00000226600 ENSG00000125356 ENSG00000203989 ENSG00000226929 ENSG00000228517 ENSG00000177485 ENSG00000226023 ENSG00000236371 ENSG00000278646 ENSG00000158290",
"ENSG00000176774 ENSG00000176746 ENSG00000232030 ENSG00000188408",
"ENSG00000198205 ENSG00000215174 ENSG00000165591 ENSG00000198455 ENSG00000186787 ENSG00000204271 ENSG00000147059",
"ENSG00000029993 ENSG00000166049 ENSG00000183862 ENSG00000102181 ENSG00000013619 ENSG00000130032 ENSG00000171100 ENSG00000160131 ENSG00000063601 ENSG00000147378 ENSG00000102195",
"ENSG00000147099 ENSG00000067177 ENSG00000184388 ENSG00000225396 ENSG00000268994 ENSG00000198034 ENSG00000125931 ENSG00000269502 ENSG00000184911 ENSG00000275520",
"ENSG00000102081 ENSG00000176988", "ENSG00000050030 ENSG00000131263 ENSG00000187969 ENSG00000147100",
"ENSG00000102174 ENSG00000184735 ENSG00000175809 ENSG00000165186",
"ENSG00000185448", "ENSG00000101977 ENSG00000101981", "ENSG00000177504 ENSG00000205642",
"ENSG00000130150 ENSG00000101958 ENSG00000181544", "ENSG00000260802 ENSG00000077279 ENSG00000077274 ENSG00000101901",
"ENSG00000062096 ENSG00000205667 ENSG00000101825", "ENSG00000169933",
"ENSG00000134597 ENSG00000102078 ENSG00000056277 ENSG00000147262",
"ENSG00000147256 ENSG00000171054 ENSG00000147255", "ENSG00000238210 ENSG00000169551 ENSG00000186376",
"ENSG00000181191 ENSG00000130052 ENSG00000090776", "ENSG00000102409 ENSG00000182916 ENSG00000166681 ENSG00000185222 ENSG00000180964 ENSG00000204065 ENSG00000133134",
"ENSG00000197565 ENSG00000101844 ENSG00000101842 ENSG00000101843",
"ENSG00000267978 ENSG00000010404 ENSG00000123584 ENSG00000241489 ENSG00000268738 ENSG00000171116 ENSG00000269556 ENSG00000197620 ENSG00000185247 ENSG00000156009",
"ENSG00000169906 ENSG00000047230 ENSG00000182287 ENSG00000182798 ENSG00000169249 ENSG00000126010",
"ENSG00000005889 ENSG00000130741 ENSG00000123130 ENSG00000184831 ENSG00000165182 ENSG00000130066 ENSG00000174010",
"ENSG00000204086", "ENSG00000001497 ENSG00000147065", "ENSG00000076770 ENSG00000123728",
"ENSG00000169084 ENSG00000214717", "ENSG00000073464", "ENSG00000132446",
"ENSG00000126970", "ENSG00000147246", "ENSG00000173674 ENSG00000184368 ENSG00000173681",
"ENSG00000046647 ENSG00000198759 ENSG00000046653 ENSG00000196459 ENSG00000176896 ENSG00000123595 ENSG00000046651",
"ENSG00000123561 ENSG00000123572", "ENSG00000183304 ENSG00000177138",
"ENSG00000234469", "ENSG00000196368 ENSG00000224109", "ENSG00000004848",
"ENSG00000174460 ENSG00000131724 ENSG00000175556", "ENSG00000131725",
"ENSG00000087842 ENSG00000165197", "ENSG00000182220 ENSG00000234636 ENSG00000185753 ENSG00000238205 ENSG00000180182",
"ENSG00000251192 ENSG00000147121 ENSG00000065923 ENSG00000102218 ENSG00000147119",
"ENSG00000165175", "ENSG00000158301 ENSG00000198932 ENSG00000125962",
"ENSG00000196767", "ENSG00000188158", "ENSG00000198798 ENSG00000214107 ENSG00000099399 ENSG00000120289 ENSG00000169297",
"ENSG00000156697 ENSG00000122121 ENSG00000188706 ENSG00000122122 ENSG00000085185",
"ENSG00000102362 ENSG00000101811 ENSG00000174225 ENSG00000188917 ENSG00000182489 ENSG00000007952",
"ENSG00000131831 ENSG00000177324", "ENSG00000102172 ENSG00000230797",
"ENSG00000215301 ENSG00000188937", "ENSG00000147223 ENSG00000089682 ENSG00000133131 ENSG00000133138 ENSG00000165376",
"ENSG00000002586", "ENSG00000215113 ENSG00000147162 ENSG00000186810 ENSG00000215115 ENSG00000147174",
"ENSG00000186288 ENSG00000186462", "ENSG00000101888 ENSG00000176076",
"ENSG00000102128 ENSG00000133169", "ENSG00000157514 ENSG00000170935 ENSG00000147224 ENSG00000147234",
"ENSG00000165194", "ENSG00000079482", "ENSG00000147003 ENSG00000130234 ENSG00000102010 ENSG00000165195",
"ENSG00000146950 ENSG00000101850", "ENSG00000183837 ENSG00000235961 ENSG00000198681",
"ENSG00000158164 ENSG00000198908", "ENSG00000008056 ENSG00000126756 ENSG00000126759 ENSG00000078061 ENSG00000126767 ENSG00000102265",
"ENSG00000198814", "ENSG00000268009 ENSG00000165584 ENSG00000269791",
"ENSG00000173698", "ENSG00000187268 ENSG00000205542 ENSG00000123594",
"ENSG00000197976 ENSG00000196433", "ENSG00000102230 ENSG00000067992",
"ENSG00000005302", "ENSG00000125363", "ENSG00000189037", "ENSG00000269405",
"ENSG00000147394 ENSG00000147400 ENSG00000147383", "ENSG00000184788 ENSG00000124429 ENSG00000147180",
"ENSG00000072506 ENSG00000158423 ENSG00000072501", "ENSG00000102226 ENSG00000102225 ENSG00000182872 ENSG00000130985 ENSG00000130988 ENSG00000147123",
"ENSG00000122145", "ENSG00000147117 ENSG00000147124", "ENSG00000008086",
"ENSG00000187808 ENSG00000198918", "ENSG00000101940 ENSG00000068354 ENSG00000015285 ENSG00000101945 ENSG00000102317 ENSG00000102312 ENSG00000147155",
"ENSG00000130054", "ENSG00000273696 ENSG00000269586 ENSG00000270946 ENSG00000271449 ENSG00000181433 ENSG00000278085 ENSG00000169446 ENSG00000278289",
"ENSG00000196664 ENSG00000101916", "ENSG00000091482 ENSG00000185915",
"ENSG00000102158 ENSG00000226784 ENSG00000131174", "ENSG00000101846",
"ENSG00000171659 ENSG00000171657", "ENSG00000260548 ENSG00000133124",
"ENSG00000147130 ENSG00000147166 ENSG00000169562 ENSG00000147140",
"ENSG00000126012", "ENSG00000130119", "ENSG00000187516 ENSG00000165169",
"ENSG00000102109 ENSG00000126768 ENSG00000102100 ENSG00000102145 ENSG00000094631 ENSG00000068308 ENSG00000187682 ENSG00000102103 ENSG00000102096",
"ENSG00000120509 ENSG00000120500 ENSG00000147127 ENSG00000186912 ENSG00000204195",
"ENSG00000006756 ENSG00000157399 ENSG00000056998", "ENSG00000184083"
), domain = c(1L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L,
13L, 14L, 15L, 17L, 18L, 20L, 21L, 28L, 30L, 31L, 32L, 33L, 34L,
35L, 37L, 38L, 40L, 42L, 43L, 44L, 45L, 47L, 48L, 49L, 51L, 53L,
54L, 56L, 57L, 59L, 60L, 61L, 62L, 63L, 64L, 65L, 66L, 67L, 68L,
69L, 70L, 71L, 72L, 73L, 75L, 77L, 78L, 79L, 81L, 82L, 83L, 84L,
85L, 89L, 90L, 91L, 96L, 99L, 102L, 105L, 108L, 109L, 111L, 113L,
116L, 117L, 119L, 122L, 125L, 129L, 130L, 131L, 132L, 133L, 136L,
142L, 143L, 144L, 145L, 150L, 161L, 165L, 168L, 169L, 174L, 175L,
183L, 184L, 187L, 190L, 191L, 194L)), row.names = c(NA, 103L), class = "data.frame")

I changed
target <- c("1", "5", "6", "32")
selected_domains = filter(gene_domains, domain %in% target)

into
target <- over10$domain
selected_domains = filter(gene_domains, domain %in% target)

and it's the same thing that @Lynn has suggested.
Thank you very much.

You just need to learn how to do a basic subset in R.
table1 <- structure(list(n = c(29L, 11L, 10L, 10L), domain = c(1L, 5L,
6L, 32L)), row.names = c(NA, -4L), class = c("tbl_df", "tbl",
"data.frame"))

selected_domains <- filter(gene_domains, domain %in% table1$domain)

Thank you very much.


text:
Extract reads from bam files by their @RG
How could I extract reads from bam files by their read groups @RG ?, I've got one file containing all reads for 5 samples, the SM is appears NONE, So, I want to extract each sample reads by the @RG tags.

Hi and welcome, please have a think about targeting your question better. What are the genome sizes and machine? What did you try command wise?

Okey thank u for answering.

So, I use ion torrent to sequence 3 bactrial whole genome and 2 samples of fungi. I used 5 codebars but when I launched the plan I Forgot to check if I put the names of each samples. When I've got the ubam file I found only one file containing all the reads. So I Want To know if I could get each sample by its @RG or ID. Knowing that the SM tag appears "None".

Typically I use samtools for operations like this.  Specifically I use samtools view with either -r or -R flag depending on the use case.
-r STR
Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

-R FILE
Output alignments in read groups listed in FILE [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

This sounds more like a job for samtools split if you want to split out all the read groups into separate bams in one go.
http://www.htslib.org/doc/samtools-split.html
samtools split [options] merged.sam|merged.bam|merged.cram

Splits a file by read group.

Options:

-u FILE1
Put reads with no RG tag or an unrecognised RG tag into FILE1

-u FILE1:FILE2
As above, but assigns an RG tag as given in the header of FILE2

-f STRING
Output filename format string (see below) ["%*_%#.%."]

-v
Verbose output

Alternatively if you do want to proceed one read group at a time, in addition to samtools view, you can also do this with GATK using the PrintReads tool with the ReadGroupReadFilter as documented here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041415172-ReadGroupReadFilter
The command line would look like this:
gatk PrintReads \
   -I input.bam \
   --read-filter ReadGroupReadFilter \
   --keep-read-group <RG:readgroup> \
   -O output.bam

There's also a filter called ReadGroupBlackListReadFilter that allows you to exclude a list of readgroups.
https://gatk.broadinstitute.org/hc/en-us/articles/360041851151-ReadGroupBlackListReadFilter
The PrintReads tool doc is here:
https://gatk.broadinstitute.org/hc/en-us/articles/360041849911-PrintReads


text:
storing SNPs in the genotype file for each chromosome in a separate file
I have already asked this question in another forum but have not got an appropriate answer, so wondering if anyone here can help me?
I have a big file with over 3 million columns (SNPs genotypes) and 800 rows (one row for each individual).
The first row is the SNPs_ID name with the SNP position on each chromosome (see the example below please), in total there are 17 chromosomes. I want to extract SNPs for each chromosome and store them in a separate file.
I have done it with "R", but I am looking for a faster solution and I am wondering if anyone can help me to figure out another solution with bash, perl or python codes?
here is an example of my input file (SNPs on chromosome 1,2,3,4...,17)
"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"....... "Ha412HOChr17:303402"
0 1 0 0 ......0
0 1 0 0 ......0
0 1 0 0 ......0
0 2 0 0 ......0
0 1 1 1 ......1
0 2 0 0 ......0

and this is an example of my desired output for example for chromosome 1:
out.chrom1
"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"
0 1 0 0
0 1 0 0
0 1 0 0
0 2 0 0
0 1 1 1
0 2 0 0
0 0 0 0
0 2 0 0
0 1 2 2

Hi, good question but perhaps post your R code?

Are the leading spaces part of your input file? Does the second row, for example really start with 4 spaces? If not, please [edit] and remove them to avoid confusion. Also, are those really the IDs? Do they all start with Ha412HO, then Chr and number(s), followed by : and the position? Or can the Ha412HO change?

I just edited my question according to your comment. And yes, the IDs start with Ha412HOChr followed by chromosome, the : and SNP position

I believe a pandas approach would be fast enough (probably faster than any other Python approach without pandas).
After converting your table to a pandas data frame, you can do something like df.filter(regex=(".*Chr01")) for the first chromosome then write the resulting data to file with the to_csv() method.
import pandas as pd

data = pd.read_csv('input.txt', delimiter = "\t", index_col = False)

for i in ["Chr01", "Chr02", "Chr03"]:
    regex_pattern = ".*" + i
    print(data.filter(regex=regex_pattern))


text:
RPGC normalisation creates artefacts at centromere
I am studying ChIP-Seq data in HeLa cells and I've started using the RPGC normalisation of deepTool's bamCoverage. MACS2 also uses this normalisation for its peak calling.
I am seeing a large number of peaks (like, over 50% of my peaks) at the centromere of each chromosome, and the RPGC normalised BW files also show this increased signal at the centromere.
Does this mean that my protein really is at the centromere of each chromosome, or is it some normalisation bias I'm not aware of?
Here's the general aspect of my peak calling (in green). Notice the large increase at the centromere of each chromosome.

And here's a zoom on the centromere of chromosome 6. In blue at the top is raw BW, then there's a lane for input normalised BW but not RPGC, then green is the peak calling by MACS2, and black is RPGC normalised BW.

Can you please tell me the exact formula for RPGC normalization?

@Upasana https://deeptools.readthedocs.io/en/develop/content/tools/bamCoverage.html

UPDATE
Here is some other protein ChIPs that I ran through the same normalisation:

As you can see, they don't seem to be enriched at the centromere. So unless anyone brings other contributions, I will accept this as answer.
RPGC normalisation doesn't appear to create artefacts at centromere. The peaks observed here are as real as can be.


text:
How I swap rownames in these data?
I have a microarray gene expression matrix like this by weird gene IDs in rows.
> head(mat[1:10,1:5])

> dim(mat)
[1] 39302    76
>

There is matched gene symbol for each of these identifiers in another matrix
>

> dim(matched)
[1] 23107     1
>

How I can have matched gene symbol with probe identifiers in the row names of my expression matrix please? The problem is, for one gene symbol we may
I tried
> merged <- merge(mat, matched)
Error: cannot allocate vector of size 6.8 Gb

EDITED

tmp = paste(matched[rownames(mat)],rownames(mat),sep="_")
> rownames(mat) = tmp

EDITED
This is my matrix after matching prob identifiers to gene symbol
> head(array[,1:10,1:5])
               GSM482796 GSM482797 GSM482798 GSM482799
1         OR2T6    0.0171   -0.1100   -0.0394   -0.0141
2          EBF1    0.1890    0.0222    0.0832    0.0459
3 DKFZp686D0972    1.9400    0.2530    0.3770    0.8310
4        ATP8B4   -0.1490    0.0690   -0.0637   -0.0527
5      NOTCH2NL    0.1540   -0.3880    0.2160   -0.0812
6        SPIRE1    0.2920    0.1690    0.5500    0.1430

but now for some genes I have several probes. For example for gene A I have several matched probes. So I have repetition for gene A. How I can take mean over the expression of repeated genes please and having an unique value?

It's not at all clear what format your data are in. Could you add dput(head(matched))?

You need to think about whether to take the average for all probes belonging to one gene, or simply append the gene name to the probe id.
If it is simply adding a gene name, maybe try this below, note I am assuming that all your rownames of mat can be found in matched
# check whether all rownames are in matched
table(rownames(mat) %in% rownames(matched))
# rename the rows
tmp = paste(matched[match(rownames(mat),rownames(matched)),1],rownames(mat),sep="_")
rownames(mat) = tmp

There are actually two questions.
The first one is about the memory issue. I believe merge from data.table would solve that issue.
The second one would be aggregating or summarizing the identifiers. As @StupidWolf mentioned, you need to have a rule, summing them up, averaging them, ... For this second part one way is using dplyr:
merged_data %>%
  group_by(probe_colum) %>%
  summarise_all(mean or sum or ...) %>%
  mutate(new_column = paste(probe_column, gene_column, sep = "_"))

The same thing can be done faster with data.table:
merged_data[, lapply(.SD, mean or sum or ..., na.rm=TRUE), by = probe_column ] [, new_column := paste(probe_column, gene_column, sep = "_"),]

.SD means subset of data and your subsets are your groups based on your probes.

hey @haci, what you proposed should work. Adding one column to the original expression matrix should not lead to memory issues. I think the problem might be coming from not being able to match the probe id to rownames of the matrix


text:
How many types of DNA methylation are there?
My current understanding is:

bases A or C can become methylated due to interaction with a methyltransferase
these methylations only occur if the A or C is surrounded by a specific sequence of base pairs
different types of DNA methylations occur in different species at different rates

Is there a comprehensive list of which sequences of base pairs allow methylations to occur?
I'm currently aware of:

dcm-methylation: CCAGG or CCTGG
dam-methylation: GATC
EcoKI methylation: AACNNNNNNGTGC or GCACNNNNNNGTT
CpG methylation: CG (or CHG or CHH)

If I'm understanding you correctly, by "types of DNA methylation" you mean "nucleotide contexts where DNA methylation occurs".
This is going to be a function of the methyltransferase proteins involved in the methylation process, and this largely depends on the organism, and even cell type being studied.
Generally, what has been observed is:

The phenomenon of DNA methylation (methyl group addition to A/C nucleotides) has been observed in both eukaryotes and prokaryotes
The nucleotide context of DNA methylation differs between the organism, and even the cell type being studied, due to the methyltransferase proteins involved in the methylation process
Methylation of RNA, predominantly in N6 position As

In humans:

The DNMT1 methyltransferase is referred to as the "maintenance" methyltransferase, due to its observation of binding to hemi-methylated CpG sites (methylated C on one strand and not the other)
The DNMT3A/B methyltransferases are referred to as the "de novo" methyltransferases for their affinity of binding to naked DNA
These proteins almost exclusively methylate CpG dinucleotides in human cells
There are cases of non-negligible methylation of cytosines in CHG and CHH nucleotide contexts (brain tissues, for example)

The proteins and their effects on DNA methylation are similar for mice, as well.
Plants tend to have similar nucleotide contexts of methylation, but have much more CHG and CHH methylation than mammals.
In prokaryotes, again, the nucleotide context is dependent on the proteins.
Like you mentioned, we have that:

The Dam methyltransferase methylates the N6 position of the A in GATC
The Dcm methyltransferase methylates the C5 position of the second C in CCAGG and CCTGG
The EcoKI methyltransferase methylates As in AACNNNNNNGTGC and GCACNNNNNNGTT

NEB has a good list of its restriciton enzymes, and how DNA methylation affects their binding/cutting, which also gives some info on the type of methylation in cells.
As for a "comprehensive list", you'll likely have to find some relevant review papers in your organism(s) of interest and compile on yourself.
But also think about how you're phrasing your question.

Is there a comprehensive list of which sequences of base pairs allow methylations to occur?

I don't think that many papers studying DNA methylation think of the nucleotide context allowing methylation to occur so much as observing the context in which methylation occurs and how it relates to proteins in those organisms.

Yes, I was referring to nucleotide context where methylation occurs. Thanks for the great response and explaining differences between organisms/cells, but my original question remains unanswered (i.e. in approximately how many contexts have methylations been found across species?). Five to ten, or possibly hundreds? I seem to only see the above 3-5 ever mentioned. Is that because this is still an active area of research? Are these discovered through experimentation with the above restrictive enzymes?


text:
Doing plot with this data
I have normalized expression values of some genes in some patients and the how many days they have survived from their diagnosis with cancer like below

EDITED

I have taken mean of the expression values for each patient and I have how many days each patient survived since diagnosis as below
>

I done so but nothing I got
> fit <- survfit(Surv(time) ~ gene, data = km)
> print(fit)
Call: survfit(formula = Surv(time) ~ gene, data = km)

                  n events median 0.95LCL 0.95UCL
gene=-0.333140816 1      1    859      NA      NA
gene=-0.307846735 1      1    347      NA      NA
gene=-0.303559694 1      1   1339      NA      NA
gene=-0.290518776 1      1     61      NA      NA

The question is if these genes increase survival time or not but likely I am doing wrong

By your help

> fit <- coxph(Surv(time) ~ gene, data = km)
> print(fit)
Call:
coxph(formula = Surv(time) ~ gene, data = km)

       coef exp(coef) se(coef)     z     p
gene 0.8664    2.3783   1.4520 0.597 0.551

Likelihood ratio test=0.35  on 1 df, p=0.5543
n= 69, number of events= 69

For plotting I guess I do need another information because
 library("survminer")

> ggsurvplot(fit, data = km)
Error in ggsurvplot(fit, data = a) : object 'ggsurv' not found

EDITED

I  divided patients based on the median of survival days to Up and Down category like below but still I am failing to visualize that
 fit <- survfit(Surv(Time, Status) ~ Gene, data = km)

> ggsurvplot(fit, data = km)
Error in data.frame(..., check.names = FALSE) :
  arguments imply differing number of rows: 70, 0, 140

Any help please?
Thanks

In general, survival analysis can be said to be composed of two steps; Cox regression, with which you calculate the "hazard ratio" based on your variables, and a "Kaplan-Meier (KM) estimate", which is used to visuazlize the data.
Here is a nice tutorial for doing survival analysis with the survival and survminer packages. The latter includes ggplot2 kind visualization of the KM curves.
For the Cox regression, you can use continuous(*) or discrete variables, in your case these would be the normalized expression value of your gene and categories of your gene respectively. You can categorize the expression level of a gene as low vs high, or low vs medium vs high, based on different quantiles.  Another approach to categorize your expression data would be to find an "optimal" cut point with maximally selected rank statistics; from the link:

This is an outcome-oriented methods providing a value of a cutpoint
  that correspond to the most significant relation with outcome (here,
  survival)

Basically you search for the cut-point where your data is split into two such that the two resulting parts of the data show the most separation/difference in terms of survival.
* Can a Cox Proportional-Hazards Model be built only with continuous predictors?

Thank you but I have about to 130 genes here for which I want to know if they together are related to survival or not; So should I take mean of their expression?

In a scenario where you have a "gene module", yes. Basically you would be creating a new variable (mean expression of your module genes) and use this in your Cox regression. If the aim would be looking at the synergistic effects of genes to survival, that is a whole different story for which you would need to use interaction terms in your Cox regression.

Sorry I edited my question because I faced a problem

The survival object should be created with time and status/event (i.e. if the patient is dead/alive) like Surv(time, status), you lack the status variable  in your data.

Sorry because I don't know if these patients alive/dead; All I know is survival days from GSE19417 accession :(

Can you try fit <- coxph(Surv(time) ~ gene, data = km) instead of fit <- survfit(Surv(time) ~ gene, data = km)? coxph seems to be way when using continuous variables. I am adding a link to the answer so it does not get lost here int the comments.

And regarding lack of patient status in your dataset, without that info, the model assumes an event has happened for all, the following is from ?Surv: Although unusual, the event indicator can be omitted, in which case all subjects are assumed to have an event. That is why you are not getting an error but please pay attention to the consequences of this assumption on your results.

For the KM curves, you should use categorical data, so your options are using i) survfit (with categories of your variable) for Cox and KM or ii) coxph with your continuous variable for Cox and survfit with categories of your data for KM.

Sorry I have a categorical column to data for visualization in my main post but I am getting error

fit <- survfit(Surv(Time, Status) ~ Gene, data = km), here the Status is not about your variable of interest, it is whether or not the event has happened. You can try fit <- survfit(Surv(Time) ~ Status, data = km). With this, you attempt to see the effect of your gene's expression categories, high and low specifically. BTW, it might be misleading to choose Status as a variable name for the categories of your gene.


text:
Human Cell Atlas - cell annotations
Human Cell Atlas Preview Datasets have been available for a while now (there was some discussion about that earlier). However, although the raw data is available, cell labels are not (for example, which cell type each cell belongs to). That seems like the most valuable part of this resource. I've heard of people using that type of information, but I have not been able to find it on the website. Is that available somewhere?

I only have peripheral knowledge of HCA so take this with healthy scepticism, but I believe that some aspects of the analysis pipelines for the project are still being finalised and harmonised, as there have been many new developments in the last year or so. It is possible that while some of those working on the data are using interim cell labels, the final labels have not been decided on thusfar.


text:
Smith-Waterman Pairwise Local Alignment Algorithm
I have written a R code for Smith-Waterman (SW) algorithm.
I wanted to check the results of my code with online examples. In all cases -except one- I found the same alignment. So there is this one example -let's say example A*- for which I found a quite different alignment result: in the example they give three equivally good alignments whereas I found only one that corresponds to no none of those three alignments (we have the same score, though).
I am now trying to understand whether example A is wrong or if my code is not working properly.
Therefore, I have two requests:

Does anyone have a dataset of the result of SW used on two sequences?
I could then compare these results with the output of my code.
Does anyone have whatever way to compute SW algorithm? Then, I could give that person the example A, which could then try this example A, and give me the output.

This would be very helpful if someone could help, I have already spent so much time on my code, without finding anything that could explain those different outputs.
Thank you!

Example A: Sequence A = "TGCTCTGT" and Sequence B = "TACTATCAT"

There are many implementations of the Smith-Waterman algorithm. However, they are not all guaranteed to give the same output because although it's an optimal local alignment algorithm, it still depends on the chosen scoring scheme (i.e. substitution matrix and gap penalties).
With this in mind, you could try EMBOSS Water online:
https://www.ebi.ac.uk/Tools/psa/emboss_water/
Or in Biopython, the Align.PairwiseAligner() (which uses SW by default when aligner.mode is set to local):
https://biopython.org/DIST/docs/api/Bio.Align.PairwiseAligner-class.html

In R, you can use pairwise Alignment in Biostrings.
## Nucleotide global, local, and overlap alignments
library(Biostrings)
s1 <- DNAString("ACTTCACCAGCTCCCTGGCGGTAAGTTGAT")
s2 <- DNAString("GTTTCACTACTTCCTTTCGGGTAAGTAAAT")
# needleman
pairwiseAlignment(s1,s2,type="global")
# smith-waterman
pairwiseAlignment(s1,s2,type="local")

I am not very sure how to implement it from scratch. If you face more problems with this example, I can have a go at it.

are there faster methods to do this than the Biostrings implementation in R? Using R, are there hash or FFT methods?

it's not related to this question (which does smith-waterman) maybe you can try https://rdrr.io/rforge/NameNeedle/man/needles.html


text:
Confidence Interval with Wilcoxon Test in Python for log-normal Distribution
I have a couple arrays of dN/dS scores, and I would like to calculate the confidence interval for each array of data. dN/dS scores are not normally distributed but are log-normally distrbuted, so I want some non-parametric approach for finding confidence intervals.
An answer in the CrossValidated forum suggested that in R you can use a Wilcoxon test to find the median-based confidence interval, a function that I do not think the python scipy.stats.wilcoxon() package does. Is there anyway of doing something similar in Python for the confidence intervals of log-normally distributed data?
For a example short array of dN/dS scores, here is the dN/dS scores of ATPases:
array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])

I am not at sure about that, I would be amazed if SciPy couldn't do this. However in R is a base function. How many columns, comparisons? What are you measuring?

I guess you mean confidence interval of the mean (or median). I looked at the array you provided, even after log it does not look normal-distributed. Use the bootstrap for non-parametric values. I modified the code byJames Triveri, see example below, and you can use the mean or median. The textbook by Efron is a good reference for the method, which is also used in R package boot. There's also awebsite which should be ok.

import numpy as np

def bootstrap(data, n=1000, func=np.mean,p=0.95):

    sample_size = len(data)
    simulations = [func(np.random.choice(data, size=sample_size, replace=True)) for i in range(n)]
    simulations.sort()
    u_pval = (1+p)/2.
    l_pval = (1-u_pval)
    l_indx = int(np.floor(n*l_pval))
    u_indx = int(np.floor(n*u_pval))
    return(simulations[l_indx],simulations[u_indx])

x = np.array([0.02288081, 0.44170839, 0.10549733, 0.17515196, 0.09449279,
       0.07110412, 0.00893079, 0.23485109, 0.14533192, 0.05449631,
       0.10281173, 0.02113355, 0.05157087, 0.01705113, 0.02651948,
       0.09352947, 0.05828018, 0.20528157, 0.02873843, 0.0141598 ,
       0.11262881, 0.06337332, 0.13689815, 0.10557703, 0.04452507,
       0.05046484, 0.40241456, 0.05939199, 0.24423569, 0.05042892,
       0.164257  , 0.03408215, 0.04737262, 0.01037463, 0.01246448,
       0.02170375, 0.0241773 , 0.05136936, 0.02393366, 0.18913979,
       0.15781334, 0.06448557, 0.04355384, 0.02821125, 0.08015629,
       0.10985432, 0.06074574, 0.15775976, 0.05678278, 0.03749782,
       0.05518756, 0.00770479, 0.21248167, 0.08005044, 0.16307954,
       0.05783565, 0.05907416, 0.07044622, 0.13227131, 0.01627556,
       0.10859962, 0.08149819, 0.05600647, 0.16098728, 0.15183062,
       0.05202344, 0.01769589, 0.00789287, 0.07777749, 0.1324942 ,
       0.12734709, 0.17146938, 0.03890857, 0.1296019 , 0.085146  ,
       0.05602965, 0.02708089, 0.11807038, 0.07848828, 0.03291032,
       0.36302533, 0.07800343, 0.06551307, 0.04676282, 0.04765273,
       0.08060882, 0.06339636, 0.03349833, 0.01224308, 0.1481316 ,
       0.31738452, 0.15690855, 0.0693822 , 0.020425  , 0.02909208,
       0.03499225, 0.03019904, 0.13722717, 0.14403507, 0.01257245,
       0.02223452, 0.07068784, 0.0544813 , 0.08738558, 0.02884046,
       0.10549474, 0.06695546, 0.01341142, 0.09440411, 0.11840834,
       0.08558889, 0.2688645 , 0.28313546, 0.15127967, 0.01463191,
       0.1728421 ])

# for mean,95% confidence interval, do
bootstrap(x,1000,np.mean,0.95)
# for median, 95% confidence interval, do
bootstrap(x,1000,np.median,0.95)

This is a good answer and it would certainly work. This is the 95% CI of the mean and it's not the Wilcoxan. If the OP is using W it suggests the data is separated into columns and they are trying to determine the difference because different groups.

Essentially the question isn't clear enough at present

I agree with you @MichaelG. You can replace np.mean with np.median, its the same concept. Yes, let's hope VeritatemAmo provides more info.


text:
Which correlation method to compute the correlation score between different clusters of Sc-RNAseq data?
I'm analyzing single cell rna-seq data and trying to compute the correlation score between different clusters. Wondering how to choose the correlation method("pearson" (default), "kendall", or "spearman") of the Cor function?
I'm a new comer in bioinformatics and know little about these arithmetics. Can anyone give me some advice or show me the code to decide which method to choose for my data matrix?
Below is my gene expression matrix: rows are gene names,and columns are clusters.I use the default Peasrson method to do analysis,but I really don't know if it is right.

CorOb.cor.exp <- as.data.frame(cor(CorOb.av.exp))
CorOb.cor.exp$x <- rownames(CorOb.cor.exp)
neupeptide.cor.df <- tidyr::gather(data = CorOb.cor.exp, y, correlation,CorOb.cor.exp$x)
a <- ggplot(neupeptide.cor.df, aes(x, y, fill = correlation)) + geom_raster()+scale_fill_gradientn(colors = c("blue", "white", "red"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                   size = 12, hjust = 1))+
  coord_fixed()

a

Cross-posted: https://github.com/satijalab/seurat/issues/2198

I don't think github issues are the best place to ask for general help such as this.

Disclaimer: This is not a complete answer but more a comment, posting as an answer because of character restrictions for comments. Also, I am not a statistican but a molecular biologist who did quite some reading/courses on stats.
Two scRNA-seq tools, namely SingleR and RCA (Reference Component Analysis, Li et. al, 2017) compute correlations of bulk transcriptomic data to single cells; for this they use Spearman correlation and Pearson correlation respectively (The two tools use the same bulk transcriptome reference). I am using these two tools for different purposes so I did not really compare their results, however, single cell clusters that I compute using the reference components (calcuated with RCA) is in good agreement with single cells labels computed with RCA (adjusted rand index of 0.81).
In addition to that, please take a look at scMatch: a single-cell gene expression profile annotation tool using reference datasets, they claim (scMatch computes bulk - single cell correlations like SingleR and RCA):

From this comparison we found Spearman’s correlation coefficient
  outperformed Pearson’s correlation for all six libraries.

This question does not have a definite answer. Rather, the answer depends heavily on the assumptions you wish to make regarding your data.
Given a pair of random variables $x$ and $y$, Pearson's rho has a value of
1 when there is a linear relationship between $x$; i.e. $y = x\beta + c$, where $c$ is a constant, and there is no noise.
Note that the slope, $\beta$, and intercept, $c$, of this relationship
do not affect the correlation coefficient. With Spearman's rho, we first compute the rank
of each element in $x$ and $y$. We then calculate the Pearson correlation between
these ranks.
Kendall's tau correlation coefficient also deals with the ranks of $x$ and $y$,
but is slightly more complex. The most relevant difference in this case is likely
that tau is better in the case of ties, and may be more robust for low sample sizes.
In most cases, Kendall's and Spearman's are likely to give very similar results
If you wish to assess only linear relationships between expression profiles,
Pearson's is the most appropriate. If you wish to allow for non-linear effects,
Spearman or Kendall is most appropriate (Kendall if the data contains ties).

Pearson correlation evaluates two continuous variables. Spearman correlation is often used to evaluate ordinal variables. Kendall correlation is a nonparametric test for two continuous variables.


text:
read and matching pattern with python
I have read the content of a text file into pandas and needed some help matching the pattern.
Here is the pattern, where a can any number greater than zero and b can be zero or any number greater than zero:
father=a|0 || 0|a
mother=a|0 || 0|a
daughter1=b|0 || 0|b
daughter2=a|0 || 0|a
daughter3=b|0 || 0|b
son1=a|0 || 0|a
son2=b|0 || 0|b

Here is the function that reads the content of text file to pandas
import pandas as pd
def read():
    varian=pd.read_csv('file.txt',sep = "\t")

father=varian["FATHER"]
mother=varian["MOTHER"]
daughter1=varian["DAUGHTER1"]
daughter2=varian["DAUGHTER2"]
daughter3=varian["DAUGHTER3"]
son1=varian["SON1"]
son2=varian["SON2"]

The content of the file looks like this
HROM    POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO    FREQUECY
1   1226852 G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=12   0.005494505
1   1847936 C   T   0/1 1/1 1|1 1|1 1|1 1|1 0|1 AC=1;AN=2184    0.000457875
1   2428427 C   G   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AC=4;AN=2184    0.001831502
1   2515616 G   A   1/0 1/0 1/0 1|1 1|1 1/0 1/0 AC=1;AN=2184    0.000457875
1   3801895 T   C   1/0 1/0 0|0 1/0 1/0 0|0 1|1 AC=10;AN=2184   0.004578755

Is that really what your file looks like? Doesn't the first line start with #CHROM? Does it really start with HROM? Also, please look at our formatting tools to see how to format your posts correctly and avoid confusion.

Sorry for the delay, busy week, done below

This is not a bioinformatics question, but a basic pandas question, so SO is a better place for this... But your odd looking VCF file can be parsed as follows without even using regex:
import pandas as pd

def splitter(cell, key):
    parts = cell.split('||')
    if key in parts[0]:
         p = parts[0].split('|')
    else:
         p = parts[1].split('|')
    sp = [x.strip() for x in p if key not in x]
    return float(sp[0])

def read():
    splitter_a = lambda v: v.split('||')[0].split('|')[1])
    varian=pd.read_csv('file.txt',sep = "\t")
    a_splitter = lambda c: splitter(c, 'a')
    b_splitter = lambda c: splitter(c, 'b')
    varian = varian.assign(father_a=varian.FATHER.apply(a_splitter)\
              .assign(father_a=varian.FATHER.apply(b_splitter)\
              .assign(mother_a=varian.MOTHER.apply(la_splitter)\

etc.
display(varian) or print(varian.columns) will show that you have new columns called father_a etc.

pandas regex might still help, but re isn't needed. The concept might work. ... It's a pandas solution for sure. I like the style. I'd construct a mask and then filter the columns. They are wanting to exclude homozygotes

The immediate issue is the separator can be '/' as well

Its easy. You want to remove homozygotes and here is the solution for the FATHER class. Note the MOTHER class will comprise homozygotes in the answer because you are only searching for FATHER. To perform an intersection one approach is just to concatenate the mask.
import pandas as pd

genedict = {\
            'FATHER': ['1/0', '1/1', '0/0', '1/0', '0/1'], \
            'MOTHER': ['1/1', '1/1', '1/0', '1/0', '0/0']}

df = pd.DataFrame(genedict)
mask = df.FATHER.str.contains(r'1[|/]0') + df.FATHER.str.contains(r'0[|/]1')
df[mask]

OUTPUT
  FATHER MOTHER
0    1/0    1/1
3    1/0    1/0
4    0/1    0/0

Note, this is not your original dataset, its a dictionary to demonstrate the principle and note all your other columns will carry across unless you want to 'drop' them. Finally, pandas has its own regex system, regardless the traditional r'' code is still needed


text:
Is there any way to align ChIP-seq reads to telomeres?
I know that telomeres are highly repeated sequences, but is there any way to retain any reads that map to these regions (on HG38)?
I recently managed to find some protein binding to centromeres, which are also mainly repeats. Therefore I wondered if there was any way to pick up signal on the telomeres.
My team and I are investigating binding of a protein complex to DNA and the experimental team reported that they saw it a bit on telomeres, so we were wondering if there was any way to measure this through bioinformatics. If not then it's no problem, we would just like to know.

If it's hard-masked (NNNNN) in the genome, no tool can map against it in the genome. Think about instead mapping to repeat regions or repeat consensus sequences from repeat masker.

This paper https://www.biorxiv.org/content/10.1101/728519v1.full examines the mapping of some human telomeric sequences, and indicates that there is something like 130kb missing at the telomeres in the HG38 assembly.  One option might be to take the HG38 assembly, and add some missing telomere sequences as an additional pseudochromosome to the FASTA file.  Then if you used this modified reference to align your reads, you should soon get a feel for the level of telomeric reads in your data set.

As far as I understand from this question and a related previous one, you have performed a ChIP-seq analysis and saw peaks at centromere regions. In that case, you can use samtools to extract reads that are aligning to a given genomic region:
samtools view input.bam "Chr?:XX-YY" > output.bam

That's not really what I am looking for. I am looking for ways to map reads onto telomeres, and then maybe see peaks from that. Edited my question slightly for clarity.


text:
What do the symbols mean in minimap2's gap cost equation?
Heng Li's paper on minimap2 is understandable to me right up until the moment that it gives this definition of gap cost:

2.2.1 Alignment with 2-piece affine gap cost
Minimap2 performs DP-based global alignment between adjacent anchors in a chain. It uses a 2-piece affine gap cost (Gotoh, 1990):
$$\gamma_a(l) = \min \{ q + |l| \cdot e, \tilde{q} + |l| \cdot \tilde{e} \}$$

Suddenly, out of nowhere, we've got five variables used without prior introduction: $l$, $q$, $e$, $\tilde{q}$, and $\tilde{e}$.
Heng doesn't state their meanings, and based on a quick skim and text search of the cited paper by Gotoh, that paper doesn't use these variable names anywhere either.
What does this equation mean?

Here is the translation:

$l$ = gap length
$q$ = initial gap penalty for shorter gaps
$e$ = short gap extension penalty per every additional gap
$\tilde{q}$ = initial gap penalty for longer gaps
$\tilde{e}$ = long gap extension penalty per every additional gap

To understand the first of the two expressions being minimized, see the Wikipedia section on affine gap penalties.
Notice the similar notation structure
$$A + B \cdot L$$
where

$A$ is the known gap penalty,
$B$ is the gap extension penalty, and
$L$ is the length of the gap

To understand the other variables with the tildes on top, see the text below Equation 2 (emphasis added to key explanation):

On the condition that $e > \tilde{e}$,
  it applies cost
  $q + |l| \cdot e$
  to gaps shorter than
  $[(\tilde{q} - q)/(e - \tilde{e})]$
  and applies
  $\tilde{q} + |l| \cdot \tilde{e}$
  to longer gaps. This scheme helps to recover longer insertions and deletions (INDELs).

In other words, we make the assertion that we should penalize longer gaps differently than shorter gaps.
As a bonus, here's a visual understanding of why we use the minimization function on those two expressions using Python and some made-up (and probably not biologically useful) penalty terms. Note the requirement that $e > \tilde{e}$.
# Import libraries
import numpy as np
import matplotlib.pyplot as plt

# Setup
l = np.array(range(0, 15, 1))
y_short = eval('10 + 5*l') # Q=10, E=5
y_long = eval('20 + 3*l')  # Q=20, E=3

# Plot
plt.plot(l, y_short)
plt.plot(l, y_long)
plt.legend(["Short Gap Penalty", "Long Gap Penalty"])

Using our example and the expression $[(\tilde{q} - q)/(e - \tilde{e})]$, it will apply the first expression in the minimization function (the one without the tildes) when the gap length is shorter than 5.
$$
[(\tilde{q} - q)/(e - \tilde{e})]
=
[20 - 10)/(5 - 3)]
=
10 / 2
= 5
$$
We can see from our plot that the two lines intersect at 5, which we've verified through some arithmetic.


text:
Filtering fasta files by ID
I have this code to filter some IDs from a fasta file:
for line in file_open:
    line=line.rstrip()
    if '16S' in line:
        header=line
        header=header[0] + header[1:]
        seq=file_open.readline()
        seq=seq.strip()
        print(header + '\n'  + seq)

however, this code doesn't get me all content of each ID sequence. I get this result:
>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC

When I should get this:
>ID1
AGTTTTGCCGAATTCCTTTAGGATTTTGTGAGGAACATTGTGTTAAATGTATTTATAGTTCTACTAATTC
TGCGAAGGTAGCATAATCATTAGCCTTTTGATTTGAGGCTGGAATGAATGGTTTGACGAGAGATGGTCTG
GGGCCGTGGTATACTGACCATGCGAAGGTAGCATAATCATTAGCCTCTTGATTTGAGGCTGGAATGAATG
GGGCCGTGGTATAC

Thanks

Does this answer your question? Remove/delete sequences by ID from multifasta

Well, the code you wrote prints the seq on a single line (and without any newline afterwards), are you sure that's not the only problem?

@KonradRudolph I think the question is more about writing the fasta file well formated.

@KamilSJaron No, the question isn’t about formatting. Check the output in the question again. It isn’t all on one line, it’s truncated (as in, only the first line is read and written, subsequent lines are omitted).

I thought that the truncated line could be due to copypasting. Now I see, it's not about the formatting of output, but expected formatting of the input.

You're assuming 2-line FASTA format, but the sequence can (and typically does) span multiple lines. You can use Biopython to build a simple parser:
from Bio import SeqIO

for record in SeqIO.parse("file_name.fa", "fasta"):
    if "16S" in record.description:
        print(record.format("fasta"))

This was really helpful, thank you so much


text:
correlation between transposon and gene
I have two matrix identical in N° of rows and columns. In the first matrix rows indicate transposable elements and columns  indicate samples. the matrix show the degree of expression of each element in samples:
M1 <-

            elements          ova       testes       retina   optic_lobe       suckers
1   LINES/3           1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   SINES             1.836210e+02   355.808323   577.807528   341.048786   99.902702
3  DNA TRANSP         9.181051e+00   178.661200  4387.148414    53.907711   353.413028

In Second matrix rows indicate genes and columns indicate samples. the matrix show the degree of expression of each gene in samples:
 M2 <-
            elements          ova       testes       retina   optic_lobe  suckers
1   OCBIM_22028242mg 1.473559e+03   962.953589   673.903449   558.879946   947.631992
2   OCBIM_22008718mg 1.836210e+02   355.808323   577.807528   341.048786   99.902702
3   OCBIM_22009482mg 9.181051e+00   178.661200  4387.148414    53.907711   353.413028

I have to make an element-gene correlation. I have to check, for example, if the element A expressed much in the tissue X has the corresponding gene A expressed very, little, or nothing in the tissue X
fOR example I'm expecting this result:
M1,2 <-
                elements          ova       testes       retina   optic_lobe  suckers
    1   TE/ gene1                  0.3          1          0.5         -1         0.3
    2   TE/gene2                   0.1        0.6            0         -1         0.9
    3   TE/gene3                     1          1          0.8        -0.2       -0.3

Hi @jonny, you can only find the correlation between M1 and M2. In other words, you see how well each TE in each tissue predict each gene. you do cor(M1,M2)

Maybe you can rephrase the question if you mean to do something else?

ok, I'll rephrase it. for each row, I would like to find the correlation transposable element - gene in all tissues

Ok I think I roughly get what you want to get. Below I simulate something that's like your data because you provided two identical matrices.
#3 TE, 10 tissue, random values
# number of tissue
M1 = matrix(rnorm(30),3,10)
rownames(M1) = paste("TE",1:3,sep="_")
colnames(M1) = paste("Tissue",letters[1:10],sep="_")
## make gene matrix, , random values
M2 = matrix(rnorm(30),3,10)
rownames(M2) = paste("Gene",1:3,sep="_")
colnames(M2) = paste("Tissue",letters[1:10],sep="_")
## here we make TE1 highly correlated to Gene1
M2[1,] = M1[1,] + rnorm(10,0,0.3)
## here we make TE2 highly anti-correlated to Gene2
M2[2,] = -1*M1[2,] + rnorm(10,0,0.3)
#here we calculate for each TE, its correlation with each gene.
#credit to @haci
COR = apply(t(M1),t(M2))
# you can look at the correlation
COR

What you need is the last three rows. Remember, you can only calculate correlation between the gene and the TE, and you do this across all tissue. For example, if gene1 and TE1 has a high correlation, this tells you TE1 expression can predict gene1 expression in most tissue.
It's meaningless to calculate correlation between one data point and another. You need more than 1 value. For example, if you do cor(10,5), it doesn't make sense.
Hope this is what you are trying to do..

cor(t(M1), t(M2)) will suffice, you don't need an apply() call as If x and y are matrices then the covariances (or correlations) between the columns of x and the columns of y are computed.

You are right. Thanks @haci!

thank you very much