text:
How to fix RepeatMasker.lib.nsq missing error when running RepeatModelor?
I have run RepeatModelor and I don't see the classified output file. I got this error.
Missing /home/mudith/.conda/envs/braker/share/RepeatMasker/Libraries/RepeatMasker.lib.nsq!
Please rerun the configure program in the RepeatModeler directory before running this script.

I think I have to run configure file. I found it. But it doesn't run when I just try
perl ./configure
Can someone please explain me how to run this? I'm new to this.
Also I found this in github. But I don't understand.
cd $(dirname $(which RepeatMasker))/../share/RepeatMasker
# ./configure downloads required databases
echo -e "\n2\n$(dirname $(which rmblastn))\n\n5\n" > tmp && ./configure < tmp

I don't understand the third line.

What that code snippet is doing:
That snippet you found on github (where is that from exactly?) looks like it tries to run RepeatMasker's configure script, feeding it some options that are needed by RepeatModeler.  The cryptic third line generates and dumps a block of text to a file named tmp:

2
path/to/bin/rmblastn

5

And then passes that in to RepeatMasker's own configure script as though it were entered on the keyboard.  It's the equivalent of hitting enter to accept a default, entering "2" to select the RMBlast option, typing in the path to rmblastn, hitting enter to set RMBlast as the default, and then typing "5" to save and quit.
That all looks fine, but I think this might be a problem with your (and mine, see below) RepeatMasker install rather than RepeatModeler itself.  These tools are new to me but I'm interesting in seeing this work too, so I gave a shot at troubleshooting it here.
The way I was able to get it working (note I assume you're working with conda packages here) was:

pushd $CONDA_PREFIX/share/RepeatMasker
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib to "automatically create the old RepeatMasker.lib FASTA style database" as the script's comments put it.  NOTE: you might want to first delete or move the old empty file rather than alter it because of how conda shares data between environments.  If you update it in place (which > will still do), you'll alter the file for all conda environments using the package, and conda may complaint about unexpected file sizes when checking packages.
perl configure to do RepeatMasker's setup including some makeblastdb commands using RepeatMasker.lib.  (Note this is RepeatMasker's configure script and not RepeatModeler's.) You can type in the options manually or you could probably use that code snippet you found.
popd and run RepeatModeler/RepeatClassifier and hopefully it finds the BLAST library files now

All the gritty detail:
I see a conda path in your logs.  Are you using the bioconda packages named repeatmodeler and repeatmasker?  It looks to me like the repeatmasker package in particular includes an empty RepeatMasker.lib file and that's the cause of the trouble.  I don't see any way that RepeatMasker.lib is set up automatically via the configure script, and the BLAST database files like RepeatMasker.lib.nsq that RepeatModeler expects are built from that (empty) file.
The Please rerun the configure program in the RepeatModeler directory message is a red herring, I think, in this case because the $CONDA_PREFIX/share/RepeatMasker/Libraries/RepeatMasker.lib file does exist, just not the other files for a database it assumes should exist.  If I instead run configure for RepeatMasker and update its script to not hide error messages, I see this coming from makeblastdb:
Building a new DB, current time: 10/26/2020 16:10:07
New DB name:   /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
New DB title:  /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib
Sequence type: Nucleotide
Keep MBits: T
Maximum file size: 1000000000B
BLAST options error: File /home/jesse/miniconda3/envs/example-repeatmodeler/share/RepeatMasker/Libraries/RepeatMasker.lib is empty

And it's true, that's an empty file provided by the conda package (the metadata in $CONDA_PREFIX/conda-meta/repeatmasker-4.1.0-pl526_0.json shows size_in_bytes: 0 for that one so that's apparently "correct" for the package.)
I saw a reference to that filename in a buildRMLibFromEMBL.pl script in the RepeatMasker share directory, and running that produces a FASTA file:
./util/buildRMLibFromEMBL.pl Libraries/RepeatMaskerLib.embl > Libraries/RepeatMasker.lib

Then ./configure sees a non-empty .lib file and builds the database, and RepeatClassifier finally sees the resulting database files.
But why is this needed?
Shouldn't RepeatMasker come with a non-empty FASTA file if it's needed by its own configure script?  I can see bioconda's build.sh doesn't do anything in particular to handle that file, but, maybe it should?  This is probably where somebody who knows these tools better can chime in.)

Hi Jesse, thank you very much for this great explanation. Yes, I'm using the bioconda packages named repeatmodeler and repeatmasker and you are correct. The repeatmasker package includes an empty RepeatMasker.lib file. I found the buildRMLibFromEMBL.pl script and I ran it as you said. Finally got the RepeatMasker.lib file and by running the configure file I got the RepeatMasker.lib.nsq file. Thank you so much for helping me. Now I'm gonna run RepeatModelor again. I hope it'll work this time and I'll get the final classified file. Thank you again.

Glad that helped!  As of five days ago the conversation about the bioconda packages is ongoing, so the process might change in the near future.  I'd just keep an eye on that github issue to see if they streamline it more.


text:
How to find adapter sequence
I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279)
I want to run cutadapt, but how would I find the adapter sequence so I can run cutadapt?

The definitive answer is the sequencing institute/lab. They know what protocol/chemistry etc. was used.
If you don't have access to that a number of tools check for known adapter sequences. Run e.g. FASTQC, which will tell you the proportions of adapter sequence found. Tools like TrimGalore can also autodetect the most common adapters.
UPDATE:
I checked 2M reads from one run of your dataset and it seems to me that it was uploaded already preprocessed:

Check all your samples, but I don't see a need to do any preprocessing if all of them look like this.

Thank you! @BastianSchiffthaler. That was helpful. I wanted to make matrix table of these read, how can I join columns of these read into one table in this case? because I tried to make them collection, it failed. also when you uploaded them did you upload them as fastqsanger.gz or as different format?

I updated the question with the photo of table that want to make (it is an example not the real data that I linked) @BastianSchiffthaler

I would say this is really a separate question, but essentially the easiest pipeline is: Download the fastq.gz, run salmon with a human reference (there are tons of tutorials for this). In R, use the tximport package to create a count matrix.

Thank you @BastianSchiffthaler. I used join column on collection in galaxy, why it did not work? I want to do it in galaxy.

Sorry, can't help you there, I don't use Galaxy. Please submit a new question for this so you get help.

source: Stack Exchange

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

sort: votes


text:
How do I create Protein Ligand images
I was wondering how I can create images showing the bonds and the residue codes involved in the bound Protein Ligand complex obtained after Molecular Docking as depicted in the image below.

Pymol, in either its commercial or open source edition, will be sufficient to make figures such as your examples. Pymol recognises "polar contacts" for the intermolecular interactions and has labels that can be used for residue codes.
It is also scriptable and is widely used.

thank you very much for getting back to me and helping me out.

Yes, PyMol is one of the best programs for rendering IMO. Usually, you would create folder with your PDB file and a file that contains all the commands to create your image like shown here. You can create more than one file. Because PyMol is based on Python you can use all python syntax here. E.g. you can create multiple images and with for loops you can create series of images that you can combine to a movie later on.
Other programs to check out would be VMD. I use VMD for quick videos and images. The PyMol images I like better most of the time.
Or newer software like Samson modeling platform. I am not sure about the licensing conditions. Have not tried it yet.
Biopython does not support the generation of 3D molecular images. It is more a library for sequence data.

Great answer. On the biopython point, NGLView is Jupyter notebook visualiser which can accept an impressive multitude of structures, including Biopython structures —however, this is totally not suggested for paper figures (as the OP wants).

@soren thank you very much for getting back to me and helping me out.

The other contributors have given great suggestions. I thought I ought to mention Michelaɴɢʟo, which is an online PDB file/PDB code/PyMOL-PSE/gene-name visualiser, where one can create an interactive page, add a description with special protein links (prolinks) that change the interactive protein view and then share the URL securely.
It has a modal to make the selection easier, but uses the NGL selection language and to show a ligand, one would add a residue type prolink selecting ligand. If the neighbourhood-showing radius setting is greater than 1.0 Å (default is 3.0 Å) one will likely see hydrogen bonds with the neighbours.
One can try it without registering (free), but cannot keep the page for longer than a month or edit the description.
Disclaimer: I am the author.

Thank you very much for getting back to me and helping me out.

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


text:
How to identify to each scaffold a read belongs to, inside a .sam file?
I have a fasta file assembly and combining it with the raw reads we produced a .bam file which I converted to .sam .
The .sam information lines look like this:
A00321:42:HLLVYDSXX:2:2302:6153:3505    99      NODE_1_length_3415511_cov_137.721502    16      60      128M    =       607     742     CGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTCCATCATAGCACCATTG        FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        RX:Z:TTGAGGGTATAGTAGT   QX:Z:FFFFFFFFFFFFFFFF   TR:Z:GACACCG    TQ:Z:FFFFFFF    BC:Z:AGTTGCAG   QT:Z:FFFFFFFF   XS:i:-10        AS:i:0  XM:Z:0  AM:Z:0  XT:i:1  RG:Z:over_1kb:LibraryNotSpecified:1:unknown_fc:0        OM:i:60

Separated by mandatory fields it would be something like this:
QNAME: A00321:42:HLLVYDSXX:1:1644:2248:3881
FLAG: 99
RNAME: NODE_1_length_3415511_cov_137.721502
POS: 1
MAPQ: 60
CIGAR: 1S127M
RNEXT: =
PNEXT: 536
TLEN: 386
SEQ:  ATCGGGTCTGACACCGCGATTAGTCCGGCCAAATCGCCGTCGAGCGCAATGAACATAACGGTCTTGCCCTCAGCGCGCAGCGCATCGGCCTTGGCGTCGATTGTGGAGTGCTCGACGCCCATGATGTC
QUAL: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

I'm actually interested in the meta data. I want to know how the RX: and BC: fields are distributed across the scaffolds in the original assembly.
I imagined the .sam file already contains the information about the assembly used to produce it. If I'm wrong, I'm sorry and please correct me, I'm just assuming.
What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read_ID,Scaffold_ID,Read_Position_Inside_Scaffold,RX,BC
Then I want to use that database to analyse the distribution of RX and BC inside each scaffold.
That's what I want.
Ultimately what I'm trying to do is evaluate the quality of my assemblies based on the Barcode distribution.
I'm good at programming and parsing, I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

What I want to do is, for each read in the .sam file, I find out its position in the assembled scaffold, and I record, Read_ID,Scaffold_ID,Read_Position_Inside_Scaffold,RX,BC

I'm just having trouble figuring out, where, inside the .sam file, can I find the scaffold and scaffold position of each read.

You already listed all those:

ReadID -> QNAME
ScaffoldID -> RNAME
Position -> Pos

Assembly evaluation is a wide topic, but other software you can use includes BUSCO, QUAST or something like LTRi. Here's a more in-depth guide by the SciLife lab in Sweden. Usually, when aligning reads back to an assembly, you care about the mapped insert size distribution compared to the theoretical insert size distribution, which tells you about collapsing or exapnding regions (i.e. structural misassemblies). This is something that tools like FRC_Align do.


text:
Loss of predictive power of polygenic risk score when dataset contains missing variants
I am trying to calculate polygenic risk scores (PRS) scores for a new dataset. This dataset does not have all the variants that the PRS score needs. The PRS score I am interested in has 40 variants, but in my dataset only 20 of those variants are present. What is the correct way of assessing the decreased predictive power of the PRS score due to the missing variants?
I assume it is not as easy as summing the effective weights of those variants I have and comparing to the full sum (all variants needed for the PRS).
Thanks a lot and sorry for newbie question!

In a sample with all 40, determine the accuracy of the PRS. Then in the same sample, remove the 20 you are missing in the other set, and then determine the accuracy with the missing variants. This will give you an idea of the level of attenuation. However, you may not have access to such a dataset. In that case, I would recommend imputing your data before calculating a polygenic score. (It's likely that the GWAS the PRS is from used imputation itself.)

Which software package around you using, PRSice? Generally missing values comprising half the total observations are extremely difficult to calculate using regression analysis (or many statistical analyses for that matter.
One way forward is to include (an) example(s) of all 40 variants into your data set and perform the analysis with and without the added data, then compare this against a data set where all 40 variants are observed and your samples are excluded. You might find the missing values have low probability of PRS and this would then explain your results.
Missing analysis in regression is a specialist area of statistics and is not for the faint hearted, so a work-around is the best way forward.


text:
How can I add for several bed files the header : track type=narrowPeak name=“narrowPeak” preferably in python ,can handle with R
I want to create custom tracks from these files
I can add the line : track type=narrowPeak name=“narrowPeak” manually by opening it with text editor:
track type=narrowPeak name=“narrowPeak”

but I wonder if there is no simple way to do it with code- to open all files that end with bed and add them this line.
Didnt find this option in bed tools, If it is possible to add the name extracting it from the file name it would be amazing.
Many thanks

I don't think you need to go to R or Python to do that. In bash
set -euxo pipefail
for f in *.bed
do
  bn="$(basename ${f})"
  on="${bn%%.bed}_with_header.bed"
  echo 'track type=narrowPeak name="narrowPeak"' > "${on}"
  cat "${f}" >> "${on}"
done

EDIT: I just read you want to set the name based on the filename.
EDIT2: And also to add an arbitrary directory
set -euxo pipefail

INDIR="/path/to/bedfiles"

for f in $(find "${INDIR}" -name "*.bed")
do
  bn=$(basename "${f}")
  on="${bn%%.bed}_with_header.bed"
  dn=$(dirname "${f}")
  echo "track type=narrowPeak name=\"${bn}\"" > "${dn}"/"${on}"
  cat "${f}" >> "${dn}"/"${on}"
done

Thanks, never worked in bash before. How can I set the origin directory? In addition, can you please if it's not too complicated to add a short comment for each line so it will be easier for me to google it?

You can change the working directory of your terminal with cd <dir>. You can also add the directory to the for-loop for f in <path_to_dir>/*.bed

I added this line : INDIR=/MacHD/Users/walkerlab/Documents/ChIPSeq”, kept everything else the same and it claims that the folder doesn't  exist            ++ find $'/MacHD/Users/walkerlab/Documents/ChIPSeq?\200\235' -name '*.bed'
find: /MacHD/Users/walkerlab/Documents/ChIPSeq”: No such file or directory

It looks like there is some byte sequence in your path: ChIPSeq?\200\235. Can you check for non-ASCII characters there?


text:
Are codons in RNA layered? Are we misinterpreting RNA codons?
I am analyzing nucleotide base-pair patterns in RNA and DNA, and had a thought about RNA and DNA (Let me first state though, I am not a biologist; I am an algorithmatician, employing a sort of recursive sequential neural network).
In any event, I was wondering if we are misinterpreting the way that we look at codons in RNA.  Are we?  Are we missing something when we interpret them the way that we do?  A picture is worth a thousand words, so I made a powerpoint slide of my question and screen captured it.  It highlights the specific question I am asking, which is: are codons layered?
Why is this important to me?  When I look at patterns, I know for a fact that unless it is the end of the sequence (a stop codon) that if I have a U-nucleotide, then the next two nucleotides cannot be AA, AG, nor GA.  Additionally, I know that unless there is a U-nucleotide, the sequence will definitely not end in two or less nucleotides.  This is what got me to wondering if there indeed is layered information in the combination of nucleotides.  Please see the picture.
So if we have (as in the picture) an RNA sequence:
AUGGUCAGUCCAUAA
We might interpret that as 5 codons {Methionine, Valine, Serine, Proline, STOP}
But what if there are actually 13 codons: {Methionine, Tryptophan, Glycine, Valine, Serine, Glutamine, Serine, Valine, Serine, Proline, Histidine, Isoleucine, STOP}

Read up on "reading frames" - that concept addresses the scenario you bring up here. Every coding sequence has 3 frames: 123 456 789, xx1 234 567 89x and x12 345 678 9xx. Usually, the one producing the longest product is chosen as the "correct" reading frame.

They know it's this way through experimentation. They can identify proteins and they have found proteins corresponding to coding using 3 non-overlapped letters and not the other way. It's not interpretation but observation.

Thank you Ram (that's really insightful), and thank you juanjo because that was another question I had about where the 3 comes from.  Your response is helpful because one of the difficulties in researching it beyond assigning 1's and 0's is language that helps you research prior patterns.

This relates to the concept of the reading frame, which is the point at which you start dividing 3-letter sequences into codons. As you point out, there are multiple ways to divide a string of amino acids into triplets - in fact, there are six unique ways, found by shifting the start letter by one three times, and doing the same but reading in reverse.
Typically, there is one commonly used frame which can be transcribed and turned into protein, called the open reading frame. But it is not unheard of for a single piece of DNA to encode multiple, overlapping genes, depending on where you start reading it. So, the sequence you show cannot encode 13 amino acids in a single string, but it can encode three different strings, one of which has 5 amino acids and the other two which have 4 amino acids each.
See Multiple Reading Frames at Wikipedia, which shows the following example of a single piece of human DNA which encodes two different genes, MT-ATP6 and MT-ATP8.

This is an awesome answer!  Thank you.

To add complexity to this discussion, there are examples in some mitochondrial genomes where the stop point of a gene is based on the polyA-attached transcript, rather than the genomic sequence. I worked on a paper about Nippostrongylus brasiliensis mtDNA, and needed polyA completion to explain the transcript expression I observed:
https://doi.org/10.12688/f1000research.10545.1
There are plenty of other examples of gene transcription happening in a situation other than one-exon-per-fragment: transcription that's on both sides of the DNA, exon sharing for different genes present on the same strand (they're usually in the same frame, but there's the MT-ATP6/MT-ATP8 example that @Nuclear_Hoagie pointed out), micropeptide expression from smaller chunks of genes. Some also have some degree of stop-codon read-through, so the stop codons are not a strict stop in all cases. Some viruses are also particularly good at exploiting transcription / translation complexity.

Thank you for opening me up to a whole new world of information.  This is all very interesting, and there is a lot to learn.


text:
Splitting data table up based on separate file of gene names
I am very new so bare with me. I have a file called data (.csv that I turned into a data table) with 2 columns (geneid, normalized counts) and about 9000 rows. I have a second file called Xgenes of genes that I am interested in pulling out of the first table (1 column- geneid, about 50 rows).
Is there a way to make 2 new files- I want one with all the geneid/normalized counts minus the ones in the second file (so minus about 50 genes). Then I want to make another file that only has the geneid/normalized counts of the 50 genes in the second file.
Currently I am doing this manually which does indeed work but if there is a way to do it through R it would make my future analysis much easier.
Thanks!

What have you tried? You can use df1$col1 %in% df2$col1 and !(df1$col1 %in% df2$col1).

I would look up awk or forward your R code for loading a dataframe, which will give you access to the above code. Once you've a bit of code to share thats when help is more available.

Come on, what happened to "friendly to New Contributors?"

@Megan As you are new please note that it is actually good practice to provide some reproducible examples. You can use dput() on your data to get a copy/paste ready version that you can share here. That makes it easier to understand how your data look and therefore easier to come up with solutions.

There are plenty of efficient and smarter solutions, but if you have so small tables, you can just use the very basics of R to do what you need:
# you need to adjust the following two lines to fit to your file formating and file names
genes_filt_table <- read.table(..., col.names = c('geneid'), stringsAsFactors = F)
norm_counts <- read.table(..., col.names = c('geneid', 'norm_counts'), sep = ',', stringsAsFactors = F)

genes_to_filt <- norm_counts$geneid %in% genes_filt_table$geneid
filt_genes <- norm_counts[genes_to_filt, ]
non_filt_genes <- norm_counts[!genes_to_filt, ]

write.table(filt_genes, 'filtered_genes.csv', quote = F, row.names = F, sep = ',')
write.table(non_filt_genes, 'non-filtered_genes.csv', quote = F, row.names = F, sep = ',')

This being said, there are like milions of more efficient solutions out there. I provided one that is using as basic R functions as possible, as it seemed that you would like to see a simple R solution.


text:
What are the coordinates of the extended major histocompatibility complex (xMHC) in GRCh38?
Is there an accepted boundary of the extended major histocompatibility complex (xMHC) in human? Specifically, I am interested in the boundary coordinates for Genome Reference Consortium Human Build 38 (GRCh38 / hg38).
I'd like to use these coordinates to query MyGene.info to retrieve a list of genes in the xMHC. For example, King et al 2019 used the following definition:

we excluded drugs with nonhuman or xMHC gene targets. We considered xMHC to include HIST1H2AA and KIFC1 and all genes between them (Chromosome 6 25.7 Mb-33.4 Mb).

But I was curious whether there was an authoritative definition of the xMHC region in humans that would specify down the the base pair rather than just megabase pair.

This is an extremely good question, these genomic coordinates of course exist, the only thing is MHC ... would include mice as well as humans.  The HLA.is very well defined but  I think you need to ask an immunologist because it is common knowledge, but I don't think any one here knows.

Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13
which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.
The "extended MHC" defined here [1] is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The xMHC is defined in [1] as the regions between HIST1H2AA and RPL12P1 genes, which would give you chr6:25726063-33400644 in Grch38. This is very similar to your HIST1H2AA - KIFC1 definition (chr6:25726063-33410226). Of course if you're trying to reproduce a result, probably better to use the original paper's definition.
[1] https://www.nature.com/articles/nrg1489

Thanks! Sounds like the 2004 Nat Rev Genetics article titled Gene map of the extended human MHC that you link to is the most authoritative source for the xMHC boundaries. I see in Figure 1, the xMHC goes from HIST1H2AA (extended class I subregion) to RPL12P1 (extended class II subregion). The flanking regions in grey are not part of the xMHC.
Note that RPL12P1 is contained within KIFC1, such that KIFC1 starts before and ends after RPL12P1 in coordinate space.

Here is a MyGene.info query that returns all NCBI gene IDs and symbols in between chr6:25,726,063-33,410,226 for GRCh38. I like these boundaries since KIFC1 extends past RPL12P1. For some reason RPL12P1 is not included in the MyGene results. Not sure why. It is a pseudogene, but other pseudogenes are included. The MyGene query currently returns 318 genes.

The approved symbol for  ncbigene:221613 is now H2AC1 rather than HIST1H2AA. KIFC1 is ncbigene:3833. RPL12P1 is ncbigene:729727.


text:
Mining documents for nucleotide and peptide sequences
The text mining literature has an emphasis on identifying and normalizing gene names, mutations, pathways, concepts, and so on. I haven't been able to find much, however, on methods for extracting nucleotide and peptide sequences from documents. I did find this tidbit from Aerts et al. (2008), emphasis mine:

Text was split into words and words greater than 10 characters in length with greater than 40% of characters from the capitalized DNA alphabet [ACGT] were extracted using regular expressions to isolate putative DNA sequences.
All putative DNA sequences extracted from each paper were concatenated in the order they appeared in the text into a single fasta sequence and labeled with the corresponding PMID.
Concatenation of sequences was performed to merge sequences split by line breaks in the text conversion, and because we reasoned that inappropriate joins would be reconciled at the genome level by local alignment procedures.
Extracted, concatenated sequences were used as queries to BLAST RepeatMasked versions of genome sequences downloaded from the UCSC genome database...

I'm struck by how simple the author's method for extracting DNA sequences is. But I guess with the limited complexity of the DNA alphabet, that makes sense. This approach wouldn't work with the much more complex protein alphabet.
Are simple regular expressions the text mining community's state of the art for finding DNA/RNA sequences in documents? What methods are used to identify protein sequences?

I am not sure there is a state of the art for that.
Some possible strategies I can think of:

Applying a dictionary to identify language words and discard these.
Using some context information like identifying some words in the same sentence than the sequence.
I guess the format in papers for sequences will not be the same as for text. Maybe you can find a pattern.

I would use an incremental strategy. Use the simplest algorithm then see what it recovers. If there is garbage, find a strategy to remove some garbage with a common pattern. Repeat until you are happy with the result.
Disclaimer: I have some experience in Natural Language Processing.


text:
Table error when combining counts columns
I have these files that I want to make table of (https://github.com/learnseq/RNAseqfiles01.git), the table that I want (since all the files have the same ID column),, I want merge the ID column with the count column from each file:
ID    File 1 count column     File 2 count column     File 3 count column   File 4 count column

I cleaned up the files from extra lines using:
grep -vE '^{' GSM2580324_counts.txt | grep -vE '^\\' | grep -v '^ *$' | sed '/^{/d'  | tr '\\' ' ' | sed '/}/d' | grep -vE '^__' > GSM2580324_counts04.txt

But when I try to make a table in R using this:
%R files = list.files(path=".", pattern="GSM")
%R file1 = read.table("GSM2580319_counts01.txt" , col.names=c("gene_id","GSM2580319_counts01.txt" ))
%R file2 = read.table("GSM2580320_counts02.txt" , col.names=c("gene_id","GSM2580320_counts02.txt" ))
%R file3 = read.table("GSM2580323_counts03.txt", col.names=c("gene_id","GSM2580323_counts03.txt" ))
%R file4 = read.table("GSM2580324_counts04.txt", col.names=c("gene_id","GSM2580324_counts04.txt"))
%R out.file = merge (file1, file2, by=c("gene_id"))
%R out.file2 <- merge(out.file, file3, by=c(\"gene_id\"))
%R out.file3 <- merge(out.file2, file4, by=c(\"gene_id\"))
%R outtable01 <- write.table(out.file, file = "out.file3",sep="\t", row.names = FALSE)
%R outtable01

I get this error:
RParsingError                             Traceback (most recent call last)
<ipython-input-179-2c471cf65113> in <module>
      7 get_ipython().run_line_magic('R', 'file4 = read.table("GSM2580324_counts04.txt", col.names=c("gene_id","GSM2580324_counts04.txt"))')
      8 get_ipython().run_line_magic('R', 'out.file = merge (file1, file2, by=c("gene_id"))')
----> 9 get_ipython().run_line_magic('R', 'out.file2 <- merge(out.file, file3, by=c(\\"gene_id\\"))')
     10 get_ipython().run_line_magic('R', 'out.file3 <- merge(out.file2, file4, by=c(\\"gene_id\\"))')
     11 get_ipython().run_line_magic('R', 'outtable01 <- write.table(out.file, file = "out.file3",sep="\\t", row.names = FALSE)')

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth)
   2315                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2316             with self.builtin_trap:
-> 2317                 result = fn(*args, **kwargs)
   2318             return result
   2319

<decorator-gen-160> in R(self, line, cell, local_ns)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
    185     # but it's overkill for just that one bit of state.
    186     def magic_deco(arg):
--> 187         call = lambda f, *a, **k: f(*a, **k)
    188
    189         if callable(arg):

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/ipython/rmagic.py in R(self, line, cell, local_ns)
    729             if line_mode:
    730                 for line in code.split(';'):
--> 731                     text_result, result, visible = self.eval(line)
    732                     text_output += text_result
    733                 if text_result:

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/ipython/rmagic.py in eval(self, code)
    262             try:
    263                 # Need the newline in case the last line in code is a comment.
--> 264                 value, visible = ro.r("withVisible({%s\n})" % code)
    265             except (ri.embedded.RRuntimeError, ValueError) as exception:
    266                 # Otherwise next return seems to have copy of error.

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/robjects/__init__.py in __call__(self, string)
    385
    386     def __call__(self, string):
--> 387         p = rinterface.parse(string)
    388         res = self.eval(p)
    389         return conversion.rpy2py(res)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface_lib/conversion.py in _(*args, **kwargs)
     38 def _cdata_res_to_rinterface(function):
     39     def _(*args, **kwargs):
---> 40         cdata = function(*args, **kwargs)
     41         # TODO: test cdata is of the expected CType
     42         return _cdata_to_rinterface(cdata)

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface.py in parse(text, num)
     44         raise TypeError('text must be a string.')
     45     robj = StrSexpVector([text])
---> 46     return _rinterface._parse(robj.__sexp__._cdata, num)
     47
     48

/opt/applications/python/3.6.3/gnu/lib/python3.6/site-packages/rpy2/rinterface_lib/_rinterface_capi.py in _parse(cdata, num)
    556         if status[0] != openrlib.rlib.PARSE_OK:
    557             raise RParsingError('Parsing status not OK',
--> 558                                 PARSING_STATUS(status[0]))
    559     return res

RParsingError: Parsing status not OK

What should I do?

Regarding your grep, why not just grep for what you want, i.e. ^ENSG.*? (These last few lines that htseq-count outputs - which is the software that produces your files - are really important though. You should absolutely look into them to see if something strange is going on with one or more samples.)
Your error comes from the escaped quotes \" you use here:
out.file2 <- merge(out.file, file3, by=c(\"gene_id\"))

Before, you used:
out.file = merge (file1, file2, by=c("gene_id"))

which is the correct way to quote here. On a sidenote, you don't need c() if it's only a vector with one element. You can just write by = "gene_id".
Here's a different way you could read the files:
files <- list.files(path=".", pattern="GSM")
count_files <- lapply(files, read.table, col.names = c("GeneID", "Count"))
count_matrix <- do.call(cbind, lapply(count_files, "[[", "Count"))
dimnames(count_matrix) <- list(count_files[[1]]["GeneID"], files)


text:
Number of edges in Graphical Analysis using Jupyter with Go Programming language
I'm trying to learn graphical analysis in a biological setting.
The online course suggested using the Go Programming language, so I'm new to it and help to understand why the number of edges that I read from the following code is different than from the site I'm downloading it from.
import (
    "fmt"
    "log"

    "github.com/kortschak/graphprac"
)

infile := "YeastL.dot"
g, err := graphprac.NewGraph(infile)
if err != nil {
    log.Fatalf("failed to read graph: %v", err)
}
{
    fmt.Printf("G has %d nodes and %d edges.\n", g.Nodes().Len(), g.Edges().Len())
}

The output I'm getting is: G has 2361 nodes and 6646 edges.
But from the website http://vlado.fmf.uni-lj.si/pub/networks/data/bio/Yeast/Yeast.htm it's suppose to read 7182 edges.
Thanks!

the dot file in the Git repo has 6646 edges, the pajek file in your link has 7182.
Chances are the latter has been updated since.

Thats a good point @Rohit you need to contact the authors directly specifically whether an update has occurred


text:
Why use "robust" FPKMs?
Both DESeq2 and edgeR have an FPKM/RPKM function that by default uses normalized library sizes ("robust" option in DESeq2). FPKMs have their own issues, but I thought the main benefit was to have comparable units across independent experiments. If you are adding an experiment-specific normalization factor, then why even use FPKMs?

FPKM are inherently experiment specific and can not be used to compare across samples. Let's consider the following two sequencing runs. Let $E1$ and $E2$ be the true, underlying expression in two samples of genes 1-6. Let $S1$ and $S2$ be the observed expression in our sequencing.
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 10 & 100 & 20 \\
  G2 & 100 & 10 & 100 & 20 \\
  G3 & 100 & 10 & 100 & 20 \\
  G4 & 100 & 10 & 0 & 0 \\
  G5 & 100 & 10 & 0 & 0 \\
  G6 & 100 & 10 & 0 & 0 \\
\end{matrix}
$$
Our totals are 60 counts for $S1$ and 60 for $S2$. We sequenced both libraries to the same depth, but as there are fewer genes expressed in $E2$, we just naturally capture the genes that are expressed more often. For simplicity's sake, let's assume the genes all have the same length of 1, so we ignore the "K" in FPKM. Let's also forget the scaling by a million, since that's just to make nice numbers:
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 0.167 & 100 & 0.333 \\
  G2 & 100 & 0.167 & 100 & 0.333 \\
  G3 & 100 & 0.167 & 100 & 0.333 \\
  G4 & 100 & 0.167 & 0 & 0 \\
  G5 & 100 & 0.167 & 0 & 0 \\
  G6 & 100 & 0.167 & 0 & 0 \\
\end{matrix}
$$
Now let's calculate size factors. In DESeq2, we create a pseudo-sample which is the geometric mean of all counts of each gene, so $\sqrt{10 \cdot 20} = 14.14$ for genes 1-3 and genes 3-6 are ignored as they have a member with 0-valued counts. The size factor for each library is the median ratio for each library compared to that pseudo-sample:
$$
med([\frac{10}{14.14}, \frac{10}{14.14}, \frac{10}{14.14}) = 0.71
$$
$$
med([\frac{20}{14.14}, \frac{20}{14.14}, \frac{20}{14.14}]) = 1.41
$$
So where do we end up if we divide by the size factors?
$$
\begin{matrix}
Gene & E1 & S1 & E2 & S2 \\
  G1 & 100 & 141 & 100 & 142 \\
  G2 & 100 & 141 & 100 & 142 \\
  G3 & 100 & 141 & 100 & 142 \\
  G4 & 100 & 141 & 0 & 0 \\
  G5 & 100 & 141 & 0 & 0 \\
  G6 & 100 & 141 & 0 & 0 \\
\end{matrix}
$$
So in this simple example, size factors produce a better estimate than FPKM or TPM to compare values across experiments, hence the "robustness". Please note that this is still not great for any statistical testing, and in general should be used with caution anywhere.
Finally, I'll leave some more reading material concerning FPKMs and cross-sample comparisons here:

[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o&t=32m
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package: https://support.bioconductor.org/p/80733/#80736
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839

Great explanation. Those would be normalized counts and I get why they make sense. But this does not really answer the question of why mix normalization and FPKMs.

Right, maybe the point I was trying to make wasn't exactly clear. The "M" in FPKM is attempting to correct for library size, but simply taking the column sums is too naive of an approach. Using size factors as a library size estimate is a more robust way to calculate that (hence the name of the option). It's not really mixed, it just replaces what's used as the library size term in the formula, since the default is inadequate in many respects.


text:
How to use EPIC DNAm files (.idat) to infer copy number alteration
I have access to EPIC DNA methylation profiles from approximately 40 samples. I want to use these .idat files in order to predict copy number alterations within specific promoter regions of interest. Can you recommend any R packages that could perform this task?
Through research I have stumbled upon epicopy (https://github.com/sean-cho/Epicopy/blob/master/vignettes/Epicopy.Rmd), but I'm uncertain whether this is the best package for this problem.
Cheers,
GB

This is a good question and very interesting area of investigation. The package you sourced looks okay from my 'outside the field' perspective, but I understand your concerns.

I am not too familiar with this, but conumee (https://bioconductor.org/packages/release/bioc/vignettes/conumee/inst/doc/conumee.html) seems promising


text:
Where can I find Single Cell Data with Location "Coordinates"?
Does single cell data typically have the following meta-data: the "coordinates" (e.g. on a tissue, adjacent tissues) saying where each cell in the sample was located relative to other cells? If not, is it possible to reconstruct this with other meta-data on the cells?
With the ultimate goal of working hands on with such location-tagged data, I am hoping for suggestions on the correct terms to search for, and even references to previous studies and/or easy to use public datasets.
For visual reference of the idea I have in mind, consider this picture:

I expect expression of any given gene in the skin cells in Condition 1 to be very different according to the location of the cell. For condition 2, I wouldn't expect the expression of a gene in the skin cells to be overly different, if at all. I want to see if I can formalize this idea using data mining techniques that were created for other purposes, but first I need the proper data.
edit: Here is an example of the format I eventually hope to work with.
Fake table 1: Each row represents measurements taken for a unique cell, while the last three columns have coordinates for its location.
\begin{array}{r|lllllllll}
  \hline
 cell & phenotype1 & phenotype2 & \dots & gene1.expr & gene2.expr & \dots & x.loc.coord & y.loc.coord & z.loc.coord \\
  \hline
1 & 0.69 & 1.34 & \dots & 1.91 & 0.21 & \dots & 1.12 & 0.05 & 1.09 \\    2 & 0.34 & 0.92 & \dots & 1.74 & 2.03 & \dots & 0.57 & 0.46 & 0.24 \\
   \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\
  n & 1.97 & 1.3 & 0.96 & 0.19 & 0.66 & \dots & 0.25 & 0.02 & 1.27 \\
   \hline
\end{array}
Alternatively, if the (x,y,z) coordinates of Fake Table 1 are unavailable or not reconstructable, are there datasets which help us construct an adjacency list of cell pairs according to their location, such as Fake Table 2?
Fake Table 2: An adjacency table with pairs of cells which were next to each other when the measurement was taken.
\begin{array}{r|cc}
  \hline
pair & cell1 & cell2 \\
  \hline
1 & 1 & 2 \\
  2 & 1 & 3 \\
  3 & 1 & 4 \\
  4 & 2 & 3 \\
  5 & 2 & 4 \\
  \vdots & \vdots & \vdots \\
  m-1 & n & n-2 \\
  m & n & n-1 \\
   \hline
\end{array}
I come from a statistics background with a basic understanding of this type of data from talks I've attended, but not hands on experience. I mostly just want to explore this type of data to inform my future work. I am open to any type of scRNA-seq data with cell location "coordinates," if that's the correct thing to ask.

Hmmmm could you please explain what the end goal of your research is? You seem to be requesting RNAseq of skin cancer. Would thiis be HPV associated?

Hi, thank you for the question. I come from a statistics background and want to see if certain data mining techniques I have in mind are worth developing for this type of data, and overall for my exploration of this area for future work. So I am open to HPV and any other type of data with cell location "coordinates," if that's the correct thing to ask.

Its just the images you have shown could be a HPV associated skin cancer. Okay I guess you are looking for classification based on morphology and the shift in expression would help assess the classification. Its just the question isn't really clear and skin cancers are known for their irregular morphology

I agree with you after I reread my question. I gave it an edit to be more clear. Thank you for the suggestion to look into the HPV morphology datasets!

Does single cell data typically have the following meta-data: the
"coordinates" (e.g. on a tissue, adjacent tissues) saying where each
cell in the sample was located relative to other cells?

No. Typical scRNA-seq is just capturing random cells in a tube with no additional information. The technology you are looking for is spatial transcriptomics where you are measuring RNA levels at a particular location, but that is still not on the single-cell level.

If not, is it possible to reconstruct this with other meta-data on the
cells?

There have been efforts to reconstruct spatial relationship. For example: Satija et al. 2015. There should be more recent approaches that I am not aware of.

Thank you for the clarification and for the reference! I'll read more about spatial transcriptomics.

Is any of this what you want?
https://support.10xgenomics.com/spatial-gene-expression/datasets

Yes, I believe that is what I need. I was hoping this type of spatial data would have a clearer delineation of what cell each gene expression comes from, but that seems like a hard problem the technology hasn't quite solved yet based on what I have learned.

I would suggest searching with terms such as "spatial transcriptomic" or "spatial proteomics". You could also search for platform-specific data such as "CyCIF", "CODEX", "MIBI", "slideseq", "visium" etc..
Here is a website that has curated some spatial data (http://www.spatialomics.org/SpatialDB/). We also recently released about 3TB of spatial imaging data that might be helpful. Unlike transcriptomic data that is often submitted in GEO (https://www.ncbi.nlm.nih.gov/geo/) there is still no one place where all of the spatial data is deposited and so you will find them only as supplemental data in published manuscripts. Hope that helps.

A good solution is to use Human papilloma virus (HPV) PCR (see below) or particularly RNA expression of protein E5, E6 and E7 of the viral genome in skin cells. High expression of E5, E6 or E7 would be diagnostic of skin cancer because of over deregulation of latency of the virus. This is because the viral control regions (L1/L2? - I forget) get disrupted mostly by genome integration leading to E5-E7 disrupting the human cell cycle control.
Over expression would be diagnositic because alot (most? all?) skin cancer is a result of HPV. HPV types 16 and 18 are the key, but most work done in this area focuses on cervical cancer because cervical HPV infection has a high mortality rate. I suspect skin cancer has a comparatively low mortality rate (but is of considerable imporant of course).
Alot of work has been done on the tandem arrays of HPV genomic insertion into the human geneome disrupting the regulatory genes leading to over production of E5-E7 and again is considered important in malignancy. However, its the protein expression of (E5), E6 and E7 that is critical. It is not clear whether human genome integration resulting in a viral tandem arrays that disrupt viral proteins L1/L2 will automatically trigger cancer, or whether it is merely a high risk factor (i.e. tandem array human genome insertions could still remain benign).

Thank you for your insight, this seems like a promising direction. Here's hoping I can find a PCR dataset for HPV with subsamples from many neighboring tissue locations.


text:
Discrepany between AutoDock Vina and AutoDock Tools
I have been calculating the binding Affinity between a Ligand and a Protein using both AutoDock Vina and AutoDock tools but the binding affinity estimates obtained from the 2 molecular docking softwares are very dissimilar. For instance, in the example below the binding affinity estimated by AutoDock Vina is -12.6 but for the same Ligand and Protein the Binding Affinity according to AutoDock Tools is +158.0. Why would there be such a massive difference.
Thank you

You also scored a 2.11 kcal/mol in ligand efficiency, which is too good to believe as around 2.1 is the theoretical max! The autodock vina score has the correct sign and size, whereas the ADT looks more like a Ki value. I am not sure what is happening, but a good place to start is by looking at the dlg file and have a look there —there are lots of boring bits, but if you search for "lamarckian" or "____" (which is a table) you get to the fun stuff.

@MatteoFerla I added the table from the dlg file. How do I determine the which value ADT has calculated. Is there anyway to find out if the ADT value is the Ki value. Thank you.

Your RMSD are 70 Å: this is a lot. most protein are 20-30 Å in width. -800 kcal/mol is released by TNT, so those scores are really bad. My guess is that one or two files changed. Have you touched the map files between Vina and ADT?

@MatteoFerla I haven't changed the map files. The binding affinity estimated for other Ligand Protein combinations by AutoDock Tools and AutoDock Vina are alike but for 2 Proteins i.e. Methyltransferase and Exoribonulease the binding affinty estimated by the two softwares is discrepant.


text:
Lipidomics normalization (peak area)
I have data from a lipidomic experiment comprising: 2 conditions and 5 samples in each condition; and for each lipid the value is a peak area.
I am doing a lipidomics analysis using the web app lipid ontology, and they recommend to normalize the data prior to doing the analysis,

"There are many ways to normalize data, so LION/web cannot check
whether data is properly normalized. If you don't know how to
normalize the data, we recommend normalization by expressing all lipid
species as fraction of the sum. A popular website that supports
several normalization algorithms is MetaboAnalyst."

I don't know how to normalize my data, I cant find any paper related to lipidomics normalization to know to best way to do that, and also I can't find normalization function in the MetaboAnalyst.
Thank you in advance!

Hi @Mee have you run controls common to all data sets? Are you sure this is not 'standardisation' rather than 'normalisation'?

Hi @Michael, they called Normalization all the time in the lipid ontology website. What would be the difference? Thanks

Hmmmm I'll make a when I have a moment. There is a very real difference between normalisation and standardisation from a statistics point of view. In my view this is standardisation. I need to know details about your runs. Were all runs simultaneous? same day, same time, same apparatus?

ok, I know what you mean with normalization and standardization. Still, they called it normalization in the website... The samples were not run simultaneously, but they were run in the same conditions..

Basically and briefly,

Normalisation = transforming a data set to conform to a normal distribution. Converting to a frequency is NEVER a normalised transformation
Standardisation = ensuring each point of the data set is comparable to the others, either directly or relatively depending on the calculation

The information you have received is incorrect, moreover you only convert to a frequency for a final assessment of the data and in my experience NEVER as an intermediate step. Normally data point between 0 and 1 are bad news and often require a transformation to get them away from fractions.
I don't understand the behaviour of the variance involved in metabolism and particularly lipids, however you should refer this issue to "cross-validation" Stackexchange. I did have some understanding of enzyme kinetics but in a different area.
Cross-Validation Stackexchance statisticians will understand less of the biology than I do (which is a bit worrying of course), so you will need to give them quite detailed information on the data set. They will almost certain support my concern of the advice you are given and provide a better method of standardisation. One method is to take the difference in mean and divide it by the standard deviation, you can then normalise the standardised data. They will want to know the goals of your investigation to understand who to process the data further.
What I can say is variation between runs is often an issue (in your apparatus perhaps not, I dunno) and standardising against a common control is one approach to overcome this, e.g. the control same common to all runs would be the mean about which all other data is centred.
What I suspect is that lipids are less popular than other fields, so standard protocol may not be observed.


text:
What type of Protein ID is this?
I got a FASTA file from NCBI RefSeq.
These proteins are annotated with IDs like this : protein_id=XP_019970915.1.
What type of ID is this? I would like to have them converted in TrEMBL IDs in order to do, later, a GO term enrichment but I have to know where I am starting from.
Thank you in advance

I'm not sure I understand the question, but XP_019970915.1 is an 'NCBI Reference Sequence' ID for 'erythrocyte binding antigen-181' (https://www.ncbi.nlm.nih.gov/protein/XP_019970915.1/). If you look that protein up in UniProt you can find other identifiers/info and GO annotations etc (https://www.uniprot.org/uniprot/A0A151LVZ9). Does that answer your question?

In full thanks!


text:
kraken2 to OTU table creation
I did run kraken2 and I get two kinds of output
How do I generate OTU table from these outputs for phyloseq usage?

Hi @twinkle if you don't get a response, maybe expand your question to provide details of your outputs? There are kraken users on the site.

You need to run Bracken to the Kraken2 report output to estimate abundance.

I have multiple text file as output, and text file contains data like 2 column format with thousands of rows. from these type of outputs how do I generate OTU table

you will use the --report option output from Kraken2 like the input of Bracken for an abundance quantification of your samples. E.g. kraken2 --db ${KRAKEN_DB} --report ${SAMPLE}.kreport ${SAMPLE}.fq > ${SAMPLE}.kraken where ${SAMPLE}.kreport will be your input for Bracken.


text:
I am trying to run nab to create the pdb file but it isn`t working (AMBER)
Here is the input file for nab:
molecule m;
m = arna( "gcuucuucuucuucgc" );
putpdb( "lr16.pdb", m, "-wwpdb -nocid -tr");

I am trying to run nab to create the pdb file with this command :
$AMBERHOME/bin/nab lr16.nab\./a.out

But it comes out saying :
gcc: error: lr16.nab./a.out: No such file or directory
cc failed!

I have loaded the appropriate modules , please could someone help

I have no idea about nab ... however
“$AMBERHOME/bin/nab lr14.nab./a.out”

Looking at your error lr14.nab./a.out:
It suggests you didn't put a space character between lr16.nab and a.out. If you run
$AMBERHOME/bin/nab lr14.nab ./a.out

I suspect it will work because at present it is looking for a file named lr14.nab./a.out which could be interpreted as a directory of lr14.nab. and a input file of a.out or the error suggests it is reading / as part of the file name. a.out of course doesn't exist its your specified output file name

output of ls -l
-rw-r--r-- 1 phys 89 Oct 6 12:58 lr14.nab – Peter 5 hours ago

"/nab: lr16.nab: no such file" is what it comes up with now. This is very weird as the instructions tell me to use lr16.nab. Any ideas?

what does ls -l  give you?

I will check now

-rw-r--r-- 1 phys  89 Oct  6 12:58 lr14.nab

-rw-r--r-- 1 an909 phys 145 Oct  6 12:59 README

I've modified my post accordingly. It should work now

I found out why it wasn’t working , turns out the function I was using doesn’t exist in AMBER! I’ve changed the function and it works ! Sorry for the long time it took me to reply, and thank you for your help

Try this command:
$AMBERHOME/bin/nab -v lr14.nab -o ./a.out,

then execute a.out separately using
./a.out.

You will get the pdb file of nab extension.


text:
How does cromwell and AWS batch decide which compute environment (on-demand or spot) to use?
For processing genomics workflows I make use of widdle and cromwell with an AWS back-end. I have followed all examples on the AWS cromwell tutorial. In my AWS batch console I can see there are 2 job queues (priority-cromwell and default-cromwell) and both queues can submit jobs to 2 compute environments (ondemand-cromwell and spot-cromwell). I have however no idea how I can control which jobs are put in which queue and are executed in which compute environment. I'd like to make optimal use of the spot-instances, how can I do that? Do I have to add something in the runtime block of my widdle task definitions, should I include something in the cromwell configuration file, or are for example only scattered jobs executed on spot instances?

Cromwell will send jobs to the AWS Batch queue ARN defined in the cromwell.conf file in the Cromwell server. To use Spot instances you can create an AWS Batch Compute environment configured to use Spot instances and create a queue for that compute environment. You can then set the arn of that queue in the configuration file or override the configured value using the standard mechanisms of overriding Cromwell’s configuration (workflow options, runtime parameters etc). If you change the configuration file you will need to restart the server.
Batch queues can point to more than one compute environment. Each will have a priority. Jobs in a queue will be sent to the highest priority environment until that environment has reached its max vCPU capacity, at that point jobs will be sent to the next environment.


text:
Extract sequences from MultiFASTA aligned file, by coordinates
I am trying to extract a specific sequence from a multifasta file, from each sequence in the aligned file.
The sequences look like this, and there are 32 sequences within the multiFASTA:
>NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome
ATT AAA GGT TTA TAC CTT CCC AGG TAA CAA ACC AAC CAA CTT TCG ATC TCT TGT AGA TCT GTT CTC TAA ACG AAC TTT
AAA ATC TGT GTG GCT GTC ACT CGG CTG CAT GCT TAG TGC ACT CAC GCA GTA TAA TTA ATA ACT AAT TAC TGT CGT TGA
CAG GAC ACG AGT AAC TCG TCT ATC TTC TGC AGG CTG CTT ACG GTT TCG TCC GTG TTG CAG CCG ATC ATC AGC ACA TCT
AGG TTT CGT CCG GGT GTG ACC GAA AGG TAA GAT GGA GAG CCT TGT CCC TGG TTT CAA CGA GAA AAC ACA CGT CCA ACT

Using the following code I can pull the sequence, using coordinates, from a single FASTA file, but how would one iterate this to pull the section (by coordinates) from each sequence within the multiFASTA?
from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

record = SeqIO.read("Whole genome.fas", "fasta")

with open("Spike protein.fasta", "wt") as out:
    writer = FastaWriter(out, wrap = [])
    writer.write_header()
    writer.write_record(record[21563:25384])
    writer.write_footer()

Thanks in advance!

I developed an answer as follows:
from Bio import SeqIO
from Bio.SeqIO.FastaIO import FastaWriter
import os

os.chdir("C:/Users/")

with open("Spike protein.fasta", "wt") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        out.write(str(seq_record.id) + "\n")
        out.write(str(seq_record.seq[21563:25384]) + "\n")

You can use SeqIO for the writing in addition to reading the sequences. This has the advantage that it will always write the correct fasta header, which is not always just the .id if there is whitespace, e.g. (untested)
from Bio import SeqIO

with open("Spike protein.fasta", "w") as out:
    for seq_record in SeqIO.parse("Whole genome.fas", "fasta"):
        seq_record.seq = seq_record.seq[21563:25384]
        SeqIO.write(seq_record, out, "fasta")

BTW don't put spaces in your filenames, it can lead to issues


text:
cmapPy exception
I'm trying to read gctx files with cmapPy. This is my code:
from cmapPy.pandasGEXpress.parse_gct import parse
from cmapPy.pandasGEXpress.write_gct import write
parse("C:/Users/User/Desktop/gctx/CAS9_BASELINE_LEVEL3_BING__n10721x10174.gctx")

This is the error:
(base) C:\Users\User\Desktop\hello>D:/anaconda/python.exe c:/Users/User/Desktop/gctx/readgcxtx.py
Traceback (most recent call last):
  File "c:/Users/User/Desktop/gctx/readgcxtx.py", line 3, in <module>
    parse("C:/Users/User/Desktop/gctx/CAS9_BASELINE_LEVEL3_BING__n10721x10174.gctx")
  File "D:\anaconda\lib\site-packages\cmapPy\pandasGEXpress\parse_gct.py", line 136, in parse
    num_row_metadata, num_col_metadata) = read_version_and_dims(file_path)
  File "D:\anaconda\lib\site-packages\cmapPy\pandasGEXpress\parse_gct.py", line 168, in  read_version_and_dims
    version = f.readline().strip().lstrip("#")
  File "D:\anaconda\lib\encodings\cp1252.py", line 23, in decode
    return codecs.charmap_decode(input,self.errors,decoding_table)[0]
UnicodeDecodeError: 'charmap' codec can't decode byte 0x90 in position 2440: character maps to <undefined>

I tried changing the path with double slashes or r"..." but it doesn't work. Thanks.

The error I see is,

UnicodeDecodeError: 'charmap' codec can't decode byte 0x90 in position
2440: character maps to

This suggests there is an error in the input format that the program doesn't like.
The thing that bothers me is why you have a mix of / and \ for directory paths. Windows uses \ and Linux and OS X use / . I don't Windows so I really can't comment.

It looks like you've used a character at position 2440 that the code doesn't like, this is often the case with ambiguous characters such as degenerate nucleotides or ambiguity codes for example "?" or "N". It seems to suggest its neither or these but something like ...
I'd be cautious of is the carriage returns because Windows is of course very different to Linux in this regard. \r (I think) versus \n

undefined mapping suggests is quite a specific error and therefore resolvable.

from cmapPy.pandasGEXpress.parse import parse
from cmapPy.pandasGEXpress.write import write
parse("C:/Users/User/Desktop/gctx/CAS9_BASELINE_LEVEL3_BING__n10721x10174.gctx")


text:
Is it possible to convert BAM file from one genome assembly to the other?
I Have multiple BAM files that are referenced to UCSC genome assembly GRCh37/hg19 that are read in different time frames. Now, I am planning a different studies that require assembling all the data against one assembly. Will it be possible to convert all the BAM files that I have to an upgraded version or convert the assembly to a different new assembly and how?
I am new in this and not a bioinformatician so, if you have any questions regarding my question to answer you can ask, any type of help is highly appreciated.

It is possible, with tools like CrossMap, but my recommendation is to re-align. Assemblies can change a lot between versions and your reads will undoubdedly align differently. Re-aligning is cheap, it will take a few days and you're good. Re-doing your analysis because you realize the coordinate translation messed something up will cost you much more time and energy.


text:
Multiple Alignment cost application
Is there any biological application/case where someone would be interested in the estimated total cost of the alignment between a set of sequences (genes or amino acids) without the aligned sequence itself?
Edit: to clarify by cost I mean the cost function such as the edit distance. I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such an algorithm would be useful.

UPDATE
In response to comments, it seems that the question is more along the lines of:
"what are some applications where computing summary statistics associated with a specific MSA result could be useful, even if we can't return the sequences?"
I have 2 responses to this:

in keeping with my original answer, I would say that it would still be useful to obtain some measure of identity, similarity, or pairwise distances between sequences (or any summary statistic about the relationship of the sequences). There are lots of applications for this level of data. I don't know whether this is something that can also be extracted from the computation, I'm a little unclear on what exactly the proposed computation is doing if it's getting a cost but not sequences.

Michael says below (and I agree) that an alignment cost is likely to strongly depend on particularities of the underlying alignment model as well as the data, and things get nasty when you try to account for gaps or weirdness. Nonetheless, it's possible that a cost could be used as a sort of guide heuristic for a separate follow-up alignment procedure, to indicate a portion of parameter space which is likely to yield good results. I don't know what this looks like, but it could be simple as resource provisioning for the computation (e.g. CPUs/RAM/disk). Alternately, maybe the cost result can be a little richer and actually seed some parameters in the next computation.

Again, without more detail I'm not sure what this looks like (and I'm not an alignment expert so my suggestions wouldn't be very good anyways).

The simplest application might be to come up with a quick "check" of the sequences to ask whether they are alignable in the first place. Like, as a QC step (did you forget to reverse complement?) maybe. (this assumes that the alignment in question is pretty onerous and you don't want to have to try several times.)

In summary, a lot depends on how you get that cost and what else you could get in addition to the cost. I wonder if rather than focusing on cost you can output some other richer summary of the alignment (such as e.g. sequence similarities, or gapping stats, etc.).
ORIGINAL ANSWER
The short answer is "yes".
Quite frequently we are interested in measuring "identity" or "similarity" of two sequence strings. As one totally arbitrary example, in metagenomics, quite frequently two rDNA amplicons that show 97% identity are considered to be in the same species, 95% same genus, 98% same strain, etc. So these measures are frequently useful as heuristics.
Identity itself is just a length-normalized edit distance across a sequence, after all.
But more generally, such measures of identity or similarity can themselves be understood as cost measures, in the sense that e.g. a dynamic programming alignment algorithm is trying to find an optimal set of matches between two sequences, which necessarily will maximize identity. Generally what it is actually maximizing is similarity (according to some cost matrix such as e.g. BLOSUM or Jukes-Cantor or whatever).

Thanks for the answer do you have any application in mind where the number of strings under analysis is >2. To be clearer I have some ideas for heuristics on the multiple sequences alignment problem (which is intractable) but these would not retrieve the aligned sequence, therefore, I was wondering whether such algorithm would be useful

@gcorso i see. i share Michael's skepticism below. it seems to me like it will be highly dependent on the specifics of parameterization and the actual sequences, such that it will be low generalizability. possibly it could be useful as an intermediate heuristic computation that is an input to another procedure, but I don't think it would be easy to interpret on its own.

Thank you for the answer, do you have any specific procedure in mind where it could be useful as input?

i'll just edit answer with response.

I'm not a fan personally, but it is context specific. The problem is that it is heading towards optimisation, particularly algorithmic optimisation. If the alignment is the final calculation - that is cool, but the situations are select. It depends how far you want to take it.
The rationale of my skepticism is simple - if an alignment was decided on a 'maximization' criteria, be it likelihood or some other metric - the subsequent calculations based on the alignment, such as phylogenetics Jukes Cantor, Kimura 2 and 3-parameter etc ... we enter circular mathematics. The critiera that defined the alignment is the criteria we are using to perform unrelated post hoc analysis. In truth there is a degree of circular mathematics in any algorithm used to produce an alignment, but we overcome this by optical verification precisely avoid this.
In this sense the alignment is simply to establish homologous amino acid/nucleotide positions for downstream analysis. Thus optical optimisation, whilst cautioned against in every other area of bioinformatics, is considered useful in alignments because of the risk that a logical error in the alignment is amplified in downstream analysis, potentially leading to an erroneous conclusion.
Algorithms are a problem where indels are concerned and aligning across an indel. It is very easy to make errors here and 'optical analysis' would simple introduce a blockwise deletion to avoid the uncertainty, i.e. if you are not sure just chop it out - better to lose data than risk artefacts on an equivocal judgement call.
I'm not saying don't try, I am saying don't look to overstate the applicability of the metric.

Thank you for your answer. What do you mean when you refer to a logistical error? Would it be more useful if it could produce some sort of heat map of the regions in the sequences that are conserved (again without retrieving the discrete aligned sequence)?

Sorry! I meant a logical error. An error that doesn't give a bug but could give an erroneous conclusion

Yes a heatmap is a very good idea. Highlighting areas of mutational saturation would be very useful. The theory is a bit complex, but it is one area - particularly in phylogeny  but is applicable to any analysis - where we unknowingly introduce saturated signal


text:
Extracting p-values using corrr package in R
I'm performing correlation analysis using the corrr package in R, but I'm unable to extract the p values. Does anyone know how to get the p values in corrr package in R?

If one of the solutions works for you, you can mark it as an accepted answer.

I would suggest looking at some of the documentation.
As far as I can tell, the intention of the package is to visualize and organize correlation coefficient estimates. It does not seem that the author of the package has any interest in p-value estimation.
If you are interested in computing p-values, I would suggest using a different tool, such as the builtin cor.test() function.

Here is a suggested workaround:

Generate a p-value table using one of these tools:

RcmdrMisc::rcorr.adjust()
psych::corr.test()

Then use corrr:as_cordf() to generate the tidy dataframe with the same structure as a regular corrr correlation table.
Personally, I have used the Hmisc::rcorr() package, and it works well:
Hmisc::rcorr() %>% `[[`('P') %>% corrr::as_cordf()

got it. thank you so much

I have since used this solution in my own work and I confirm that it works well.


text:
Does module size decrease as the number of samples increase in WGCNA?
We are working with a lot of samples, more than 50,000 people! We are using WGCNA to get modules and it seems like there are less number of module genes the more people we use. We are following the tutorial guidelines. We set the min modulesize = 30 and deepSplit = 4 (max sensitivity) However, we are getting very small modules with numbers almost 30 for all. We want to know if using this many samples creates an artifact and the variation between the samples may reduce the correlation seen between the genes? If that is the case, do you suggest any way we can fix that? Also do you think z-scaling expression data improves the module detection?
Thank you so much.

Weighted correlation network analysis would require a clear understanding of the a priori weightings in my opinion. The background is here

Where can access this background on priori weightings you suggested? Thanks!

The biggest difference a large number of samples makes is that you can usually decrease the soft thresholding power. If you use power 6 (or 12 for a signed network), try decreasing it to 3 or even 2 (6 or 4 for signed network). Apart from that, it could also be that the co-expression modules in your data are simply smaller than what people usually find, nothing apriori wrong with that.

Okay, thank you. I am already using 2 as a threshold. I am afraid there may be  confounding factors between the samples that are making it difficult for detection of modules. Are there confounding factors that I should be looking at removing that are known to affect module detection?

I cannot make a definite statement, it depends on the minutiae of the data. You could certainly try some noise filtering approaches although I have yet to see a bulk (not single cell) sequencing data set where these help.


text:
Extracting concatenated metadata from a column in R
I've downloaded some ChIP-seq data from www.chip-atlas.org and I want to use this to conduct some enrichment tests with genomic regions from my experiment. The .bed file generated by that website is set up so that each row represents a genomic region from a specific ChIP-seq experiment. I have the data set up as a dataframe in R.
Each row has a single column with all of the metadata, including such important information as the target transcription factor, model system, and tissue type.  These are delimited by semi-colons, but the number of pieces of information varies between rows so it's impossible to simply split by the delimiter.
The metadata is structured something like this:
name=CTCF; strain=B6; age=adult; tissue=brain
name=Mecp2; age=embryo; tissue=5t3 cells; passages=5-10

So within each row, the order and number of metadata terms varies.
I want to reshape this column into something like this:
name; strain; age; tissue; passages
CTCF; B6; adult; brain; NA
Mecp2; NA; embryo; 5t3 cells; 5-10

How can I extract and reshape this data in R? It seems like it should be possible, but the normal approaches to splitting columns by delimiters are failing me since the rows all have different numbers of items and order.
Here's a slice of the real data which is much more messy than my little example.
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

Thank you for your help!

In general you want to:

Clean your text
Split the features by semicolons
Split keys and values by '='
Create a temporary data frame with keys and value
Reshape

library(tidyverse)
library(urltools)
library(textclean)
test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%20CTCF%20brain%20ChIPSeq%3B%20Mus%20musculus%20x%20Mus%20spretus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Mouse%20adult%20hybrid%20brain;tissue=Adult%20whole%20brain;strain=BL6/spretus;chip%20antibody=CTCF,%20MillIpore%2007-729;", "ID=SRX2957352;Name=Smc1a%20(@%20AtT-20);Title=GSM2684855:%20SMC1%20ChIPseq%20in%20AtT-20%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=AtT-20%20cells;cell%20line=AtT-20%20cells;cell%20type=pituitary%20corticotroph%20cell;passages=10-20;chip%20antibody=Bethyl%20A300-055A;", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);", "ID=SRX5287318;Name=Mecp2%20(@%20Olfactory%20Nerve);Title=GSM3577864:%20GSM1827604%20MeCP2%20ChIP%20WT%20rep1%3B%20Mus%20musculus%3B%20ChIP-Seq;Cell%20group=Neural;
source_name=Olfactory%20epithelium;strain/background=C57BL/6J;genotype/variation=Wild-type,%20MeCP2(%2B/y);age=8%20weeks;cell%20type=Olfactory%20sensory%20neuron;chip%20antibody=MeCP2%20antibody%20(Diagenode,%20pAb-052-050);"))

tbl <- apply(test_df, 1, function(row){
  # Split by smicolons and clean html tags
  n <- replace_html(unlist(strsplit(row[2], ";")))
  # Decode urlencoded characters, mostly space
  n <- url_decode(n)
  # Split into keys and values
  s <- do.call(rbind, str_split(n, "="))
  # Get keys and values into vectors and clean up leading and trailing whitespace
  v <- s[, 2]
  v <- str_replace_all(v, "(^ +)|( +$)", "")
  k <- s[, 1]
  k <- str_replace_all(k, "(^ +)|( +$)", "")
  # Create final data frame with original columns (pos), and also columns from
  # keys and values
  tibble(pos = row[1], key = k, value = v)
})

# Concatenate list into tibble
(tbl <- do.call(rbind, tbl))

This tibble looks like:
# A tibble: 37 x 3
   pos     key           value
   <chr>   <chr>         <chr>
 1 3002983 ID            SRX661585
 2 3002983 Name          Ctcf (@ Brain)
 3 3002983 Title         GSM1446329: CTCF brain ChIPSeq; Mus musculus x Mus spretus; ChIP-Seq
 4 3002983 Cell group    Neural
 5 3002983 source_name   Mouse adult hybrid brain
 6 3002983 tissue        Adult whole brain
 7 3002983 strain        BL6/spretus
 8 3002983 chip antibody CTCF, MillIpore 07-729
 9 3002881 ID            SRX2957352
10 3002881 Name          Smc1a (@ AtT-20)
# … with 27 more rows

Using this, you can:
# Pivot wider
tbl %>%
  pivot_wider(names_from = key, values_from = value)

Which nets you:
# A tibble: 4 x 15
  pos   ID    Name  Title `Cell group` source_name tissue strain `chip antibody` `cell line` `cell type`
  <chr> <chr> <chr> <chr> <chr>        <chr>       <chr>  <chr>  <chr>           <chr>       <chr>
1 3002… SRX6… Ctcf… GSM1… Neural       Mouse adul… Adult… BL6/s… CTCF, MillIpor… NA          NA
2 3002… SRX2… Smc1… GSM2… Neural       AtT-20 cel… NA     NA     Bethyl A300-05… AtT-20 cel… pituitary …
3 3010… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
4 3021… SRX5… Mecp… GSM3… Neural       Olfactory … NA     NA     MeCP2 antibody… NA          Olfactory …
# … with 4 more variables: passages <chr>, `strain/background` <chr>, `genotype/variation` <chr>, age <chr>

Marvellous! This approach works perfectly, and following through the code was very helpful for my understanding of data pre-processing in general. Thank you!


text:
How do you currently track bugs in python using freely available software/websites?
I am starting a project using Python to track, archive, assign and manage bioinformatics bugs. This software will be able to be integrated into project management software or stand alone.

How do you currently track bugs/assign bioinformatics bugs using freely available bioinformatics software/websites?
What are features that you like/dislike? What are features that you would like to see in a bioinformatics bug tracker?

I am a Biology student with basic/core Python skills.

Thanks for posting this question! Unfortunately, this question is unlikely to be answered by a single best response. This type of question is probably better suited to https://www.biostars.org/

Use IPython and activate the trace in the Python debugger as follows
%pdb 1

or
%pdb on

to turn this off %pdb 0 or %pdf off
I do recommend using a formal IDE such as PyCharm or Visual Studio Code. The professional version of PyCharm is highly rated and within that use breakpoints . In Python3.7 there is supposed to be breakpoint() so you may not need an IDE, i.e. you set the breakpoint in your code using breakpoint(). In older versions of Python you can import this function via, import pdb; pdb.set_trace()
Finally check out PDB++ which is found here
Its important to note that experienced programmers often debug by simply inserting print() throughout their code and don't bother with an IDE


text:
creating a tab delimited file
I am working on a project using a fasta file. I am writing my command in nano within command-line and executing using python, also within my command-line.
I would like my command to provide me with a tab delimited file with three columns: first column should contain my sequence name, second column should provide me with my sequence length, and the third column should show the sequence itself.
I have written the following command so far within nano:
from Bio import SeqIO
import sys
for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  list = hello_fasta.split("\t")

  print hello_fasta.description
  print (len(hello_fasta.seq))

For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq
H0192X 26 FORUWOHRPPTRWFAWWEAKJNFWEJ

You can use a list and insert() to add an element in a specific order, then expand the list with *. Or you can use join().
from Bio import SeqIO
import sys

for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):
  sequences = []
  sequences.insert(0, hello_fasta.description)
  sequences.insert(1, len(hello_fasta.seq))
  sequences.insert(2, hello_fasta.seq)
  # option 1
  print(*sequences, sep='\t')
  # option 2
  print('\t'.join(map(str, sequences)))

Here's a solution using pandas if you want to save the tsv:
from Bio import SeqIO
import pandas as pd
from io import StringIO

example = """
>seq0
FQTWEEFSRAAEKLYLADPMKVRVVLKYRHVDGNLCIKVTDDLVCLVYRTDQAQDVKKIEKF
>seq1
KYRTWEEFTRAAEKLYQADPMKVRVVLKYRHCDGNLCIKVTDDVVCLLYRTDQAQDVKKIEKFHSQLMRLMELKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
>seq2
EEYQTWEEFARAAEKLYLTDPMKVRVVLKYRHCDGNLCMKVTDDAVCLQYKTDQAQDVKKVEKLHGK
>seq3
MYQVWEEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVCLQYKTDQAQDVK
>seq4
EEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVVSYEMRLFGVQKDNFALEHSLL
>seq5
SWEEFAKAAEVLYLEDPMKCRMCTKYRHVDHKLVVKLTDNHTVLKYVTDMAQDVKKIEKLTTLLMR
>seq6
FTNWEEFAKAAERLHSANPEKCRFVTKYNHTKGELVLKLTDDVVCLQYSTNQLQDVKKLEKLSSTLLRSI
>seq7
SWEEFVERSVQLFRGDPNATRYVMKYRHCEGKLVLKVTDDRECLKFKTDQAQDAKKMEKLNNIFF
>seq8
SWDEFVDRSVQLFRADPESTRYVMKYRHCDGKLVLKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
>seq9
KNWEDFEIAAENMYMANPQNCRYTMKYVHSKGHILLKMSDNVKCVQYRAENMPDLKK
>seq10
FDSWDEFVSKSVELFRNHPDTTRYVVKYRHCEGKLVLKVTDNHECLKFKTDQAQDAKKMEK
"""

# This example just happens to be a string, just load your
# fasta file using the method you're already using
example_records = SeqIO.parse( StringIO(example), 'fasta')

# Dictionary to hold the data you eventually want in the tsv
data = {"Gene name" : list(),
        "Gene length" : list(),
        "Gene seq" : list()}

# Append the necessary into the data dictionary
for record in example_records:
    data['Gene name'].append(record.description)
    data['Gene length'].append(len(record.seq))
    data['Gene seq'].append(str(record.seq))

# Convert your data into a pandas DataFrame and save as a tsv
gene_df = pd.DataFrame(data)
gene_df.to_csv("gene_info.tsv", sep = '\t', index = False)

This results in a tsv that looks like this:
$ head gene_info.tsv
Gene name       Gene length     Gene seq
seq0    62      FQTWEEFSRAAEKLYLADPMKVRVVLKYRHVDGNLCIKVTDDLVCLVYRTDQAQDVKKIEKF
seq1    106     KYRTWEEFTRAAEKLYQADPMKVRVVLKYRHCDGNLCIKVTDDVVCLLYRTDQAQDVKKIEKFHSQLMRLMELKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
seq2    67      EEYQTWEEFARAAEKLYLTDPMKVRVVLKYRHCDGNLCMKVTDDAVCLQYKTDQAQDVKKVEKLHGK
seq3    58      MYQVWEEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVCLQYKTDQAQDVK
seq4    62      EEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVVSYEMRLFGVQKDNFALEHSLL
seq5    66      SWEEFAKAAEVLYLEDPMKCRMCTKYRHVDHKLVVKLTDNHTVLKYVTDMAQDVKKIEKLTTLLMR
seq6    70      FTNWEEFAKAAERLHSANPEKCRFVTKYNHTKGELVLKLTDDVVCLQYSTNQLQDVKKLEKLSSTLLRSI
seq7    65      SWEEFVERSVQLFRGDPNATRYVMKYRHCEGKLVLKVTDDRECLKFKTDQAQDAKKMEKLNNIFF
seq8    68      SWDEFVDRSVQLFRADPESTRYVMKYRHCDGKLVLKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM

Hopefully this helps!


text:
Is there a simple command for outputting a tab delimited columns?
I am working on a fasta file and am writing my command in nano within command-line and executing using python, also within a command line.
My objective is to get my command to provide me with a tab delimited file with three columns: first column should contain my sequence name, second column should provide me with my sequence length, and the third column should show the sequence itself.
I have written the following command so far within nano:
For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq
A06842X 45 GORHQRIHQEERHEUWOHRPPTRWFAWWEAKJNFWEJQEFQEPRT

Jb4329PC 21 LDHKQAJTYOWEPLKFREEQW

Hi @KamilSJaron, I appreciate the help! The reason why I am doing this is because it makes it easier for me to easily detect which genes I could disregard if they are over  a certain length. For the command below, I tried it, but I received a "TypeError: sequence item 2: expected str instance, Seq found". I tried googling this but was not able to find anything.

Yeah, I forgot to convert the Seq class to str class before creating the line, see the edit, it should work now. Regarding the utility of such table, I guess I would always get the gene length on the fly (every time you have the sequence of a gene is really easy to get the length too in any language).

I have no idea what do you want to do with fasta files reformatted like this, but I think something like this should do the job (untested):
from Bio import SeqIO
import sys

for olig_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  name = olig_fasta.name
  seq = str(olig_fasta.seq)
  seq_len = str(len(olig_fasta.seq))

  print("\t".join([name, seq_len, seq]))

Based on OP's desired order, it should be print("\t".join([name, str(seq_len), seq]).

Yeah, you can prin Seq, but to manipulate it as a string one needs to convert it to string using str command (I edited the answer, now it should work)


text:
Is there a Python program that is able to show convergent and divergent evolution on a phylogenetic tree?
I was wondering if there was a Python program that detects convergent and divergent evolution on a tree. I am also curious if it is just better to look at a phylogenetic tree and conclude what is convergent or divergent evolution than what a program tells you.

Well in the absence of a definition of 'divergent evolution' I would simply use a definition of allopatry. If the divergence is congruent with geographic distance or geographic foci this provides some proof of divergent evolution.
In this context you are looking for something like the biogeographic BAyesian analysis such as BayArea or geographic add-ons to Beast
There isn't a singular Python library that will handle this but there is some functionality with Biopython here
from Bio import Phylo


text:
How can I subset WGS data to the level of WES variants?
I would like to compare mutational signatures1 in patients from different studies, however some studies are based on exome seq (i.e. ~20,000 coding variants) and some are from whole genome seq (i.e. ~22,000 coding variants) - is there a way to 'downsample' WGS data to better reflect WES data and effectively 'ignore' the coordinates of those ~2000 coding variants in the VCF files?

1Alexandrov, L.B., Kim, J., Haradhvala, N.J., Huang, M.N.,
Ng, A.W.T., Wu, Y., Boot, A., Covington, K.R., Gordenin, D.A.,
Bergstrom, E.N. and Islam, S.A., 2020. The repertoire of mutational
signatures in human cancer. Nature, 578(7793), pp.94-101.

Apparently this is trivial - according to a colleague you can take the WES bed file coordinates and remove any SNVs in the WGS file that fall outside those regions - but he didn't know of any software that does this in a more sophisticated manner.

bcftools would be my choice, I'm sure bedtools could do the trick, too.
Something along this
bcftools view --regions-file

or  --targets, --regions-file might require a tabix index.

I'll expand slightly on the previous answer. First print off the positions for the exome file and then use bcftools view to filter the variants from the whole genome file. You could also index the whole_genome.vcf file to make the filtering faster.
bcftools query \
    -f'%CHROM\t%POS\n' \
    exome.vcf > exome_variants.txt

bcftools view \
    -T exome_variants.txt \
    whole_genome.vcf > whole_genome.exome_positions.vcf


text:
How to remove batch effect from TCGA and GTEx data
I want to play some differentially expressed genes for gliomas cancer. I collected cancer tissue profiles from the TCGA database and normal tissue profiles from the GTEx database. I straightforward use ComBat from R package "sva", and the result shows little differentially expressed genes (DEGs) were found. How can I remove batch effects correctly?
Thanks

You can't because batch is confounded with cancer/normal. All cancer are from batch=TCGA and all normals from batch=GTEx.

Yup.  I'm very sorry, but RNASeq doesn't work like that.  It is far too sensitive to batch effects.

There is a paper that did what you are asking: Unifying cancer and normal RNA sequencing data from different sources.

This has been asked and answered here and Biostars. My take home message from many similar posts is "not to bother with correcting batch effects when performing differential expression analysis but to include 'batch' as a covariate". The latter accounts for the batch effects.
EDIT: This does not really answer OP's question, as @swbarnes2 and @ATpoint pointed out, they have a completely confounded design and the linear regression model would not accept such a model matrix.

This cannot possibly work if batch is totally confounded with sample type.

Thanks for your answer. Is it OK to ignore the batch effects?

Personally I don't think it will be OK to ignore the batch effects. If you still perform the DEA, you should be extremely cautious about the results, you will not be able to tell if the differential expression is the result of true biology or batch effect. You might want to turn to the bench to validate your results if possible.

I would at least check the documentation of TCGA and GTEx whether they at least used the same library prep kits. If not it is basically just guessing if you ask me regardless which statistics or machine/deep learning you throw at it. Then go for the genes that in your DE are really highly signiicant, and try to validate findings with independent analysis or experiments. This DE analysis here is basically just hypothesis generation, nothing more.


text:
annotation using ChIPseeker package error
I have differential binding sites object obtain from diffBind (dba.report). I am using the ChIP Seeker package to annotate the peaks but keep getting the following error: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’
I used diffloop package to add chr to the chromosome column but it did not work. Result in the same error.
Please let me know if you have any idea how to fix it. Highly appreciated
Here is my peakset
BD

  seqnames              ranges strand |      Conc Conc_condition Conc_Control      Fold      p-value         FDR
               <Rle>           <IRanges>  <Rle> | <numeric>     <numeric>    <numeric> <numeric>    <numeric>   <numeric>
    [1] NC_000003.12   75668856-75670016      * |  12.27197    12.8243533      0.76175  12.06260 1.43013e-102 1.41852e-97
    [2] NC_000024.10   11320736-11323578      * |  10.58988    11.1419124      0.76175  10.38016 1.15927e-100 5.74930e-96
    [3] NC_000004.12         37431-38634      * |   6.42233     0.0596816      8.06241  -8.00273  6.67673e-98 2.20750e-93
    [4] NC_000017.11   21678164-21686752      * |   9.26507     9.8163416      0.76175   9.05459  3.54723e-95 8.79607e-91
    [5] NC_000022.11   12692637-12694102      * |   9.64945    10.2010195      0.76175   9.43927  2.35838e-87 4.67845e-83
    ...          ...                 ...    ... .       ...           ...          ...       ...          ...         ...
  [226] NC_000004.12 131722954-131744030      * |   11.8156      12.36301     5.337173   7.02584  2.54289e-17 4.78151e-16
  [227] NC_000015.10   68134159-68136364      * |    8.2514       8.80099     0.969455   7.83153  2.28289e-16 3.73532e-15
  [228] NC_000005.10   34189082-34194027      * |   10.4671      11.01445     3.990668   7.02378  1.46745e-15 2.12829e-14
  [229] NC_000018.10       108003-108613      * |   11.2783       5.46414    12.912836  -7.44869  3.68643e-12 3.26794e-11
  [230] NC_000004.12   49709306-49710113      * |   12.0662       6.40350    13.698765  -7.29526  1.23922e-08 6.37694e-0

peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
> peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
>> preparing features information...         2020-11-12 01:08:37 PM
>> identifying nearest features...       2020-11-12 01:08:37 PM
>> calculating distance from peak to TSS...  2020-11-12 01:08:38 PM
>> assigning genomic annotation...       2020-11-12 01:08:38 PM
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’```

that was solved by changing seq level:
newnames <- paste0(c("1","2","3","4", "5","6","7", "8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","X","Y","MT"))
names(newnames) <- paste0(c("NC_000001.11","NC_000002.12","NC_000003.12","NC_000004.12", "NC_000005.10","NC_000006.12","NC_000007.14", "NC_000008.11","NC_000009.12","NC_000010.11","NC_000011.10","NC_000012.12","NC_000013.11","NC_000014.9","NC_000015.10","NC_000016.10","NC_000017.11","NC_000018.10","NC_000019.10","NC_000020.11","NC_000021.9","NC_000022.11","NC_000023.11","NC_000024.10","NC_012920.1"))

peaks_with_new_seq_level <- renameSeqlevels(peaks_with_old_seq_level,newnames)```


text:
Merging several fq.gz files or R1 and R2 classes into a single one
I have the following files:
V300066187_L4_B5RDBATtnuRAAAAA-408_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-406_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_1.fq.gz

the ending defines whether it is forward _1 or reverse _2.
I want to merge all the _1.fq.gz files into like Merge1_.fq.gz and the same for _2.fq.gz
I tried doing cat *_(number).fq.gz > merged(number).fq.gz and when I used the merges for assembly with MaSurca, it didn't detect any reads and failed. I used one pair and it worked well, but the genome size was smaller than expected, real size = 2 GB, estimated = 1.8GB.
Would appreciate how to solve this issue and use all my short reads for assembly - I need to input one file for all forward reads and one file for all reverse reads.

you can use paired reads if you define the sample names in the config file. Check the sr_config_example.txt in your current installation folder. Edit the line PE= pe 500 50 with your file names, remove the JUMP line for mate pair.

There's discussion of this topic here: https://www.biostars.org/p/81924/

Hi there, cat *_(number).fq.gz > merged(number).fq.gz should have worked. Perhaps one of the files is corrupted? You can check that using gzip -t.

What you did should work.  You don't need to decompress and recompress.
I'd check the files.  Were they trimmed?  Can you get the top, say 1000 lines of the merged files to work?  (You will need to decompress to count lines accurately)

OP seems to have two distinct V3000* prefixes - could that be the cause for the tool to fail?

Who knows?  I'm not familiar with the downstream tool, and no one can troubleshoot "it failed".  More details are needed.


text:
RAD Seq Data Analysis without barcode
I received a ddRAD Seq data from my supervisor without barcodes and restriction enzymes. I asked him for both and he said I don't need it since the data has been cleaned by the company. Now, I am kinda confused and don't know where to start.

How about having an honest face2face conversation with him rather than asking strangers in the internet. Lack of proper conversation with peers is in my experience a major source for poor production, just talk to him and ask for clarity why you don't need it and how you should proceed.

Wow Biobash, that's kinda harsh. It is cool that sometimes you unexpectedly learned something new. What Michael mentioned is pretty interested, no one has ever mentioned that or I have never heard that.
So, my approach would be that you now are doing a Poolseq study, since you can't identify individuals. The other way to try to catch the paralogues is to map your reads to a sister species. My question is what are you trying to estimate? Genetic differentation? Genetic diversity? Pop Gen structure?
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6116966/
https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.5240
https://www.g3journal.org/content/9/12/4159
One semi-robust approach to assess usefulness is thinking about your system and defining what you expect to find, let's say in terms of genetic differentiation or genetic structure, and estimating that from your data. If that hypothesis is confirmed, your dataset might be useful. You can test the same by first estimating the loci by mapping your reads to a sister reference genome, to try to detect paralogues. Since you don't have barcodes, then each "locus" is independent, and you lose power, but there is nothing else you can do.
Good luck;

The fundamental problem with RAD is that if you hit paralogue - you're done.
Identifying a paraologue isn't necessarily easy depends on the species. Essentially if you're doing pop gen you need homologues and RAD has too many situations where it can fall over without you knowing.
What happens if you've a paralogue is that you start confusing concerted evolution (pop gen of a genome) with population structure between individuals. I do agree that without barcodes identifying individuals - thats a bit tricky ;-) , however RAD an approach loaded with pitfalls.
Pop gen means being able to identify populations and on face value that will not be trivial without barcodes (I suspect they are still there), but RAD is far from an ideal tool and can give weird and wonderful metrics which have little to do with the population structure.
Its not the only problem with RAD for the purpose you are using it.

Thanks @all, I really appreciate your kind responses.

@BioBash upvotes preferred