text:
How to (efficiently) search nuccore to retrieve only annotated sequences?
I wish to efficiently search nuccore to retrieve only annotated sequences.
By "annotated" I mean that the nuccore entry contains annotation information.
An equivalent definition (if I understand correctly) is that the entry's Feature Table isn't empty. The table can be downloaded by clicking Send to in the web page of the nuccore entry.

Alternatively, the table can be downloaded using Biopython: Entrez.efetch(db='nuccore', id=nuccore_id, rettype='ft')
This definition suggests one answer to my question: Search nuccore normally, then for each entry download its Feature Table, and discard each nuccore entry whose table is empty.
The problem is that this would be very inefficient.
My best idea thus far is to add gene[Feature key] to my search term.
E.g., instead of "escherichia coli"[Primary Organism], use ("escherichia coli"[Primary Organism]) AND gene[Feature key].
This seems to work, but also feels like a quite dirty hack.
Is this a good idea? Or am I missing something crucial here?
Another idea I had was using (Biopython) Entrez.esummary(db='nuccore', id=nuccore_id), and hopefully finding there a is_sequence_annotated flag or something similarly wonderful. One of the fields Entrez.esummary returns is Flags, but I have no idea how to parse it. I asked about it in another question.


text:
Analyzing a blast result of NGS data
I am analyzing BLAST results from total RNA NGS, pair ends. This BLAST results was performed against the ssRNA virus database download from NCBI.
I filtered those read data whose paired ends uniquely matched a single virus and I have both (removing single reads).
Moreover, based on the sstart and send values from a unique virus in each library I calculated the difference (ssdif) of the lowest (sstart) of the highest (send) value that I got in the BLAST results,
for example: If I got 4 reads that match against one virus in one library (2 reads with pair ends) I got the lowest sstart and the highest send for these four hits.
Then, I calculated the coverage using the ssdif; coverage=(ssdif/genome length for that virus) * 100
After I calculate the coverage for each virus in each library based of the BLAST results I calculated the mean coverage for each genus and for each family of all viruses detected in my analysis.
Questions

Does it make sense the way that I calculated the coverage? I would like to filter my data for those viruses that have similar or higher coverage to the genus or family (previously calculated)

Based on what taxonomic level does it make more sense filter them?

Thanks in advance

source: Stack Exchange

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

sort: votes


text:
merge the matrix from computematrix of deeptools
I have several matrix from computematrix of deeptools. I need to merge two of them using "computeMatrixOperations cbind -m input1.mat.gz input2.mat.gz -o output.mat.gz"
but I am running to error "/computeMatrixOperations.py", line 521, in cbindMatrices
    hm.matrix.matrix = np.hstack((hm.matrix.matrix, np.empty(hm2.matrix.matrix.shape)))"
Is there any other way to merge the matrixes?

Please post the full error message, not just a single line.

Thanks, here is the error: https://www.dropbox.com/s/zudmn4lwo3yikw6/Capture.JPG?dl=0

You're trying to cbind two matrices of different sizes, so computeMatrix is throwing an error since it's unclear which row goes with which. You need to ensure that matrices have exactly the same number of rows and that the Nth row in each match (i.e., they should be bound together with cbind).

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


text:
awk shell script
I have a folder named quant. Inside this folder I have 27 folders named SRR8068516_quant to SRR8068543_quant. I want to edit a file inside these folders titled quant.sf and produce an output. I attempted writing a shell script to do so:
    #!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' quant.sf > equant.sf
done

When I run the script it gives me an error saying : Processing sample quant.sf8068516
awk: cannot open quant.sf (No such file or directory).
What is wrong that I am doing.
Kindly help.

You are running the command on the same file, quant.sf each time. You need to use your $fn variable instead:
#!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
    samp=`basename ${fn}`
    echo "Processing sample ${samp}"
    awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' "$fn" > equant.sf
done

However, that would also fail since it will be simply overwriting a file called equant.sf in your current directory with the results of processing each file.
I can't be sure, since you haven't actually shown us your file structure, but I suspect what you're looking for is actually this:
#!/bin/bash
for fn in /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf{8068516..8068543};
do
    samp=`basename ${fn}`
    outName="${fn/quant.sf/equant.sf}"
    echo "Processing sample ${samp}"
    awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' "$fn" > "$outName"
done

The outName="${fn/quant.sf/equant.sf}" replaces the first occurrence of quant.sf in the variable $fn with equant.sf, keeping the path and thereby ensuring that the output of each file is stored in the same directory with the name equant.sfNNNNNN. For example:
$ fn=/home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf8068516
$ echo ${fn/quant.sf/equant.sf}
/home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/equant.sf8068516

This means that your awk command will become:
awk 'BEGIN { OFS=FS="\t" } { sub("\\..*", "", $1); print }' /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/quant.sf8068516 > /home/usr/Downloads/salmon-1.2.1_linux_x86_64/quants/SRR/equant.sf8068516

Thus speaks the oracle of Awk


text:
Is there a tool that can perform a read-group-aware mpileup from a single file?
I would like to perform a samtools mpileup from a single file that contains thousands of read groups with different SM tags. I could split the bam by read group using samtools split, and then perform the mpileup, but splitting into thousands of files and then performing an mpileup from thousands of files is crazy slow on my machine.
I could imagine an implementation using pysam, but I was hoping someone has a solution already.
Bcftools mpileup already has this feature, but the feature was never implemented for the original mpileup format, which is what I need.
An example
Contents of input.sam:
@HD VN:1.6  SO:unknown
@SQ SN:chr1 LN:1000000
@RG ID:0    SM:sample_0
@RG ID:1    SM:sample_1
r0  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:0
r1  0   chr1    24  0   1M  *   0   0   G   I   RG:Z:1

When I run samtools mpileup I get this behavior:
> samtools mpileup input.sam
[mpileup] 2 samples in 1 input files
chr1    24  N   2   ^!G$^!G$    II

Whereas, what I actually want is this output:
[mpileup] 2 samples in 1 input files
chr1    24  N   1   ^!G$    I   1   ^!G$    I

For a long time I thought this was a bug in samtools mpileup, but it appears that it is a feature that was described/hinted at in the samtools documentation, but was never implemented: https://github.com/samtools/samtools/issues/599#issuecomment-604941884

been a long time since i use that.. you need the RG tag field the reads right? mpileup uses this to separate into different columns (i.e samples)

Won't GATK handle this? At the very least the old unified caller required input exactly like this.

That's a really good point! I did not realize gatk has a Pileup command. Unfortunately, that command (I tested v4.1.7.0) is also not read group aware :(

I broke down and finally scratched this itch. I have implemented a tool that performs read-group-aware text-pileup from a single SAM/BAM file. The tool is called streaming_pileupy and it's available from pypi for installation with pip. After installation the command would be:
spileup input.bam sample_names.txt

The tool is pre-alpha and missing some features. Contributions are welcome.


text:
How to modify execution orders for Snakemake?
My goal is to process two samples, S1 and S2, using two rules, step1, step2.
The Snakemake file is like:
SAMPLES = ['S1', 'S2']

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: "{sample}"
  output: "{sample}.step1.done"
  shell: "touch {output}"

rule step2:
  input: "{sample}.step1.done"
  output: "{sample}.done"
  shell: "touch {output}"

The current execution order is:
step1 on S1
step1 on S2
step2 on S1
step2 on S2

Is there a way to change the order to:
step1 on S1
step2 on S1
step1 on S2
step2 on S2

Basically, I want to complete S1, then S2.
Is there a way to specify the order of executions in Snakemake?

I think you should use the priority directive to give precedence to downstream rule(s).
Here below I give higher priority to step2 so when the first sample has completed step1 snakemake will run step2 on that sample rather than submitting another sample to step1.
SAMPLES = ['S1', 'S2', 'S3']

wildcard_constraints:
    sample= '|'.join([re.escape(x) for x in SAMPLES])

rule all:
    input: expand("{sample}.done", sample = SAMPLES)

rule step1:
    priority: 1
    output: "{sample}.step1.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

rule step2:
    priority: 10
    input: "{sample}.step1.done"
    output: "{sample}.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

(Note that without wildcard_constraints your/this code throws AmbiguousRuleException)

it seems cumbersome to create them to ensure execution order in my case

It took me a while to figure out, but it did not seem cumbersome to me.
You can re-use your rules, see below, and make Step 1 of S2 depend on Step 2 of S1 like so:
SAMPLES = ['S1', 'S2']

wildcard_constraints:
    sample = "S[12]"

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: data="{sample}"
  output: data="{sample}.step1.done"
  shell: "touch {output.data}"

rule step2:
  input: data="{sample}.step1.done"
  output: data="{sample}.done"
  shell: "touch {output.data}"

# step2 of S1 has the flag file as output

use rule step2 as step2ForS1 with:
    input: data="S1.step1.done"
    output: data="S1.done", flag=touch("flags/step2-s1-done")

# step1 of S2 needs the flag file as input

use rule step1 as step1ForS2 with:
    input: data="S2", flag="flags/step2-s1-done"
    output: data="S2.step1.done"

The problem that you are encountering, is that when snakemake determines the order to run a set of rules it first builds a Directed Acyclic Graph (DAG) of all the rules. This shows which rules are are dependant on each other, and what order they can be safely executed in, without breaking downstream dependencies.
In general, if you want tasks to execute one after another which do not share an explicit file dependency, you can use a Flag File as an ouput for rule 1 then reference it as an input in rule 2.
Unfortunately, I do not think this will help you in your case. You are expanding on the SAMPLES list and creating a wildcard. in the DAG this creates two competely separate execution branches which have no dependencies with each other. step1 and step2 are guaranteed to occur in order for each sample, however there is no guarantees between samples.
Is there a specific reason why you need all of S2 tasks to be run after S1? If so, then your best bet would be to create explicit rules for sample 1 and sample 2, using the Flag File to incidate when step2 for S1 is done.

In typical analysis, we would like to see some results for the first sample, and make sure it is correct. So I would like to see S1 processed all the way first. Thanks for offering the Flag File feature. But it seems cumbersome to create them to ensure execution order in my case.

@zhanxw your need seems to mostly center around QA of the pipeline, and making sure that its producing the right outputs. I often have to validate a pipeline before running 10s or 1000s of samples. The best way I found is to not try to gate the pipeline  like you are suggesting, but just run it once with a single sample,  then run it again with all of the remaining samples once you are satisfied.

I see your points: run one sample at a time. I can do that, but that is not convenient as I need to write another script to run snakemake on each sample. I am curious if snakemake has already some solutions. Thanks.


text:
phmm visualization
I have used HMMER search and I have created a phmm. But I would like to visualize it, not as a logo but as a phmm. I do not have Linux, so I cannot use build hmmer. Is there any web browser that I could use?Thank you!

HMMER searches databases for protein homologs. This can contribute to making sequence alignments. HMMER uses an advanced statistical approach (profile hidden Markov models).
PHMM is a profile hidden Markov models, i.e. an output of HMMER which would be used to represent a protein family.


text:
What does the number mean in an HGVSp annotation?
Let's take the example p.Arg452Pro that I got from an annotated VCF file that contains an HGVSp column.
What does 452 mean in this case ?
Theoretically, is it possible to have the same HGVSp on the same gene but on a different transcript level ?
For example:
Chromosome      Position        Ensembl_Gene        Ensembl_transcript      HGVSp
4               4613173         ENSG00052342        ENST00141456            p.Arg452Pro
4               4613176         ENSG00052342        ENST00141458            p.Arg452Pro

Thank you in advance.

the 452 amino acid of the protein?

It's the position in the protein. p.Arg452Pro means that there is a variant that changes an Arginine residue at position 452 of the protein into a Proline.
And yes, it is absolutely possible to have the same HGVS p. notation on the same gene but on different transcripts. More than possible, it is quite common. Your particular example isn't actually a real variant since that position on chromosome 4 is intronic in both hg19 and hg38, but TP53:p.R175L1 is a real example of what you describe:

As you can see in the image or by visiting the link, this variant on TP53 has the same HGVS p. notation on 6 transcripts. It also has a different p. notation on two more transcripts that are not shown in the image as well as no p. notation for a few other transcripts in which it doesn't affect the protein product.

1Disclaimer: The link goes to https://varsome.com, a free human variant search engine that is developed by the company I work for.

I actually gave a random position lol... but thanks for the explanation. Then, I have another question, what does the protein position 452 mean? Is it the number of the protein for that gene on the exons for that specific transcript ? it's not the genomic position, right ? For example, if we have a different transcript, with the position 453, can I say that the previous amino acid is 452 Arginine ?

@Law it's the position in the protein that is produced when the transcript is translated. The number refers to the protein sequence, not the genomic sequence and not the transcript sequence. The 1st amino acid is at position 1, the second at position 2 etc. Your putative variant would be at position 452, it would affect the four hundredth and fifty second amino acid. Whether that is the same for another transcript will depend on what that transcript codes for.

Do you know if clinical diagnostics are transcript-dependent, or would the same protein change on different transcripts be considered equivalent for clinical purposes?

@BradyGilg that depends entirely on the transcript and protein in question. Do both proteins have the same function? Are they expressed in the same tissue? This isn't a simple question. You might want to post a new question instead.

What terdon has answered is right. For more details you can visit the official HGVS recommendations here


text:
What is PyMOLs search path for scripts?
I know I can run my PyMOL scripts if I specify a valid path from the current working directory and also using absolute paths, but can I have a central script repository, so I don't always need to specify the full path or copy my scripts around from project to project?
Environment info:
Kernel: Linux 5.6.8-1-MANJARO x86_64
Python: 3.8.2
Python2: 2.7.18
PyMOL: 2.3.0 Open-Source

You can use pymolrc and create aliases to your scripts.
If your script is located e.g. in
/Users/tillyboy/pymolscripts/myscript1.pml

Just enter this line to your pymolrc:
alias myscript1, run /Users/tillyboy/pymolscripts/script1.pml

Now myscript1 is available for all your Pymol sessions independent of your current working directory.


text:
Is there a way to get the code from "github.mit.edu"?
everyone.
I am a postdoc who just begins to analyze single-cell RNA-seq data.
Nowadays, I found a really interesting paper (https://doi.org/10.1016/j.celrep.2018.10.047, PMID: 30404002). So, I decided to study bioinformatics by drawing the same graphs presented in the paper during the corona outbreak season. However, it was turned out that I can't get any access to the website in which code is deposited. Here is the address for the website.

Do you guys know how to get the code to learn bioinformatics and make graphs?
Thank you.

Follow the Bioconductor single-cell workflow: https://osca.bioconductor.org

It is very strange that they put their code somewhere cannot be easily accessed. If they just put it under a normal Github repo, there would not be such a problem. One way to get the code is simply to email the corresponding author. People are usually responsive if you are interested in their work.

Thank you. I just sent the e-mail to PI. I hope to get the code information.

It is not uncommon for Institutions to set-up their own private github / gitlab instances to have more control over their code and what they do with it. Furthermore, it is quite likely that this was a private repo and someone forgot to change it to public. It happens. As the saying goes, never attribute to malice what can be explained by incompetence (or absent-mindedness) :)


text:
How to query NCBI (Nucleotide database) by a feature qualifier?
I cannot seem to find a way to query nucleotide database (https://www.ncbi.nlm.nih.gov/nuccore) by a specific qualifier of the feature. For example, by /host under the feature "source".
For example, below is a part of one record, where the qualifier /host contains "Bemisia tabaci B".
How could I search for all sequences that contain the word "Bemisia" in the qualifier \host?
FEATURES             Location/Qualifiers
 source          1..351658
                 /organism="Candidatus Portiera aleyrodidarum BT-B-HRs"
                 /mol_type="genomic DNA"
                 /strain="BT-B"
                 /host="Bemisia tabaci B"
                 /db_xref="taxon:1206109"
 gene            1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
 CDS             1..1926
                 /gene="dnaK"
                 /locus_tag="C530_RS00005"
                 /old_locus_tag="C530_001"
                 /inference="COORDINATES: similar to AA
                 sequence:RefSeq:WP_014894920.1"
                 /note="Derived by automated computational analysis using
                 gene prediction method: Protein Homology."
                 /codon_start=1
                 /transl_table=11
                 /product="molecular chaperone DnaK"
                 /protein_id="WP_014894920.1"
                 /translation="MIRIIGIDLGTTNSCVAVLENGKSRVIENSEGARTTPSIIAYTD
                 DNEILIGQTAKRQAVTNPSNTLYAIKRLIGRKFNDDIVQKDIKMVPYKIVSAENGDAW
                 VQVKNRKLAPPQISAEILKKMKKTAEDYLGEKITEAVITVPAYFNDSQRQATKDAGRI
                 AGLNVKRIINEPTAAALAYGMDKNRGDKKIAVYDLGGGTFDISIIEIADVDGETQFEV
                 MATNGDTFLGGEDFDLKLINYIVKQFKLESGVDLSGDSLAMQRLKEAAEKGKIELSSS
                 QQTDINLPYITADKTGPKHLTLKITRAKLESLVEDLIKRSLLPCKKAIDDAKLTTKEI
                 DDIILVGGQTRMPLVQKKVAEFFGKEARKDVNPDEAVAMGAAIQGGVLGGDVKDVLLL
                 DVTPLTLGIETMGGIMTPLIEKNTTIPTKKTQIFSTAEDNQTAVTIHALQGERKKALQ
                 NKSLGRFDLTDIPPAPRGTPQIEVAFDLDANGILNITAKDKATGKEQSIIIKSSGGLS
                 EEEINKMIQDAETNKEEDKKFEELVQTRNQADGMIHTIRKTLKEIKNEDEKNKLESLI
                 NNLEQTLKSDDIKTIKEDLNKLTEEYSTIYQKIYSEKKPQNHNSEPAQKNAEKSSSKK
                 EGDVVDAEYEEVNNKTK"
                 ...

Best,
Milan

Doesn't it work if you just search for [host]"Bemisia tabaci B"?

@terdon Unfortunately it does not. If I try exactly as you wrote, it says "The following terms were ignored: [, ]".
If I reverse the query to be "Bemisia tabaci B"[host], then the response is "Unknown field was ignored: [host]."

Yes, it says The following terms were ignored: [, ], but the results seem to be what you're looking for, are they not?

Almost, it does include sequences that have the term in \host, but it also returns many sequences that do not have \host qualifier at all, but contain the term elsewhere, usually in the name of the sequence (DEFINITION). This is especially visible if you do [host]"Bemisia".

For this particular query, try host=Bemisia tabaci B[WORD].  However, you may also be interested in related queries like host=Bemisia tabaci B biotype[WORD].
But, I am not sure a general solution exists, as I do not believe all qualifiers are indexed.

Also, you can see how your query is interpreted on the right side under "Search details".  For the original attempt ([host]"Bemisia tabaci B"), it is interpreted as:
host[All Fields] AND "Bemisia tabaci B"[All Fields]

while the goal is to produce
host=Bemisia tabaci B biotype[WORD]

Thanks for confirming it to me, the issue is indeed that this qualifier is not indexed. In the end, I downloaded the results as GenBank file for a much wider search, and then locally imported into BioSQL database where I was able to write queries for every qualifier and get the result I needed. It takes some effort, but the flexibility of running complex queries was amazing to have.


text:
How to combine two Genome-wide Association Study (GWAS)?
I did a GWAS analysis in the past for antibiotic resistance of E. Coli and the results were interesting (matching the literature). I did a new GWAS analysis for some new samples, but the results are not in agreement with the previous GWAS and the literature. So, I was wondering if the combination of them would be in agreement with the literature, as having more samples to support an argument is always a good idea.
So, my question is how I can combine two GWAS study and what quality control steps do I need to follow?

Hi there, you got to be much more specific with your questions on this webpage. You should write the reasoning behind combining the combining and what do you want to do with the combined GWAS output. What did you consider already, what you found when you tried to research this question on your own... You can [edit] your question and add all that :-)

The first way is to combine two samples and re-run GWAS. A second way is to do a meta-analysis using GWAS results.

For combining them and re-running GWAS, how should I take into account batch-effect, as these groups belong to two different batches?
What do you mean by meta-analysis exactly?

Batches may be a problem. I do not know how to do it in E. coli, but if it's human, I would do a PCA to see if the two batches overlap. I think the batch effect is moderate in genotype data. In meta-analysis, you can just use the effect sizes and p-values from many studies to calculate a new p-value. You can read some paper in human GWAS perhaps.

What about SNP missingness for different batches? Should I remove the SNPs with high missingness for one batch and low missingness for another batch? Or Should I look at the average missingness.

Any answer, dear @Phoenix Mu?

I think if you want to merge the genotypes and re-do GWAS, you should look at overall (average) missingness. I have no idea about the SNP missingness in E. Coli. In human, it is usually a very small proportion of SNPs that are missing and should not be a concern.

I'd recommend first having a think about how you would go about combining two GWAS with the same sample set (e.g. split the samples randomly in two). Even with the same source population, you'll get different associations coming up by sampling a different sub-set of that population.
Once you can do that multiple times, you can attach a level of confidence to each SNP association within each study, which can be treated as a prior probability for the purpose of combining multiple independent studies.
It sounds like you're in a great position to be able to do this, because you've done the analysis on both groups and [presumably] have all the genotype-level data from both studies.

Dear @gringer
I am not sure if I understood you totally.
First, it is good to know that my samples in the two groups are different. Do you mean I don't stick with these two groups, rather I mix them, iteratively?!

First, stick within a single group. That way you don't need to worry about batch effects and SNP dropout. How would you compare two different GWAS taken from the same group?

But, my two groups are coming from different sources, that's why I didn't mix them from beginning.

Yes, I'm aware of that. But it's a good idea to start with something that works, and extend it to do what you want it to do. Work out first how you'd combine results from the same source. If you can't do that, there's no hope in combining results from different sources.
Once that's sorted out, see how much of that combination process can be applied to results from different sources.


text:
EGA remote download via snakemake fails due to EGA_CLIENT_ID specification
The example code provided by the Snakemake documentation for remote downloading from the EGA archive fails for me. I have defined all required environment variables, but an error comes up stating
{"error":"invalid_client","error_description":"Bad client credentials"}

I guess I am not defining the EGA_CLIENT_ID variable correctly? I am using Snakemake 5.16 and Python 3.7.2. Thanks a lot in advance for your advice.
Code:
import snakemake.remote.EGA as EGA

ega = EGA.RemoteProvider()
os.environ["EGA_USERNAME"] = "XXX"
os.environ["EGA_PASSWORD"] = "XXX"
os.environ["EGA_CLIENT_SECRET"] = "XXX"
os.environ["EGA_CLIENT_ID"] = "pyega3"
rule a:
  input:
      ega.remote("ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai")
  output:
      "data/COLO_829BL_BCGSC_IlluminaPipe.bam.bai"
  shell:
      "cp {input} {output}"

"Bad client credentials" looks more like your username/password combination was not recognized. Do they match your EGA account?

Thanks for the comment - yes the match my account, I took them straight from the CREDENTIALS_FILE which I can successfully log in with when invoking pyega3 from the command line.

Ok, that's good. Have you tried setting the environment variables before the line ega = EGA.RemoteProvider? Maybe the values are already used at this point?

Unfortunately changing the order results in the same error.

Okay, I finally solved it. The reason was the EGA_CLIENT_ID variable. In the pyega3.py script, there is a definition for the client_id for pyega3 v3.0.44 ("f20cd2d3-682a-4568-a53e-4262ef54c8f4"). After setting the EGA_CLIENT_ID to this value login was successful.
os.environ["EGA_CLIENT_ID"] = "f20cd2d3-682a-4568-a53e-4262ef54c8f4"

I think this is a bug and I filed it on Snakemake's git repository - the official documentation just requires definition of the EGA_USERNAME and EGA_PASSWORD variables, and while definition of EGA_CLIENT_SECRET will probably need to be done by the user, the client ID should be pulled automatically from the pyega3 version installed locally.
Of note, after successful login, the example download from the documentation failed with the following:
Missing input files for rule a:
ega/EGAD00001002142/COLO_829_EPleasance_TGENPipe.bam.bai

However, downloading files that I have access to worked, indicating that the paths to the files from the documentation are wrong or are access-protected.


text:
PAM50 gene expression classification
I'm looking at running the PAM50 classifier on RNA-Seq from 138 breast cancer samples. However, the R package (genefu) that's useful for this does not have a particularly helpful vignette when it comes to processing the RNAseq counts because it uses pre-processed publicly-available RNAseq data.

Starting from counts, is there a way to pre-process the data for running intrinsic.cluster.predict in R/genefu using the built-in signatures? (I probably won't have enough power to split the set into a training/test set).
The best choices for the parameters in intrinsic.cluster.predict would also be very helpful.

For differential expression (I have pre and post surgery pairs), I have used conditional quantile normalisation to get normalised log2(RPKM) values for running through R/EdgeR, so I also have these available if that is the best way to go.

EDIT: looking at the genefu vignette, this is all based on microarray data (which is essentially log2(expression)+offset), so will adding an offset to the log2(RPKM) values improve prediction for PAM50?


text:
Passing custom distance matrix to T_Coffee
I am passing a custom distance matrix to T_Coffee (written in C), as per the docs, but can't seem to get it to work.
The documentation says to pass a matrix like this using BLAST format
https://tcoffee.readthedocs.io/en/latest/tcoffee_technical_documentation.html#blast-format-recommended
# BLAST_MATRIX FORMAT
# ALPHABET=AGCT
A G C T
A 0 1 2 3
G 0 2 3 4
C 1 1 2 3
...

I've tried all combinations (with and without the leading # comments, and with and without the trailing) ..., but no matter what I keep getting this error:

Parsing_error when reading blast_matrix doolittle: A C 0

Here, my distance matrix file is titled doolittle because it is the difference between the hydropathy indices of amino acids based on the Kyte-Doolittle scale * 10 to make them all integers. I've pasted the full matrix below.
I am calling the command from linux 18.08 ubuntu terminal (although I've tried most other combinations.)
t_coffee -in=rhomboids.fasta -matrix doolittle

Any help would be appreciated. Thank you!
doolittle
# BLAST_MATRIX FORMAT
# ALPHABET=ACDEFGHIKLMNPQRSTVWY
A C D E F G H I K L M N P Q R S T V W Y
A 0 7 53 53 10 22 50 27 57 20 1 53 34 53 58 26 25 24 27 31
C 7 0 60 60 3 29 57 20 64 13 6 60 41 60 65 33 32 17 34 38
D 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
E 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
F 10 3 63 63 0 32 60 17 67 10 9 63 44 63 68 36 35 14 37 41
G 22 29 31 31 32 0 28 49 35 42 23 31 12 31 36 4 3 46 5 9
H 50 57 3 3 60 28 0 77 7 70 51 3 16 3 8 24 25 74 23 19
I 27 20 80 80 17 49 77 0 84 7 26 80 61 80 85 53 52 3 54 58
K 57 64 4 4 67 35 7 84 0 77 58 4 23 4 1 31 32 81 30 26
L 20 13 73 73 10 42 70 7 77 0 19 73 54 73 78 46 45 4 47 51
M 1 6 54 54 9 23 51 26 58 19 0 54 35 54 59 27 26 23 28 32
N 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
P 34 41 19 19 44 12 16 61 23 54 35 19 0 19 24 8 9 58 7 3
Q 53 60 0 0 63 31 3 80 4 73 54 0 19 0 5 27 28 77 26 22
R 58 65 5 5 68 36 8 85 1 78 59 5 24 5 0 32 33 82 31 27
S 26 33 27 27 36 4 24 53 31 46 27 27 8 27 32 0 1 50 1 5
T 25 32 28 28 35 3 25 52 32 45 26 28 9 28 33 1 0 49 2 6
V 24 17 77 77 14 46 74 3 81 4 23 77 58 77 82 50 49 0 51 55
W 27 34 26 26 37 5 23 54 30 47 28 26 7 26 31 1 2 51 0 4
Y 31 38 22 22 41 9 19 58 26 51 32 22 3 22 27 5 6 55 4 0

I'm sorry that I am only able to give you non specific advice, but have you tried the most recent version of t-coffee? Which version are you using? Maybe you could try a smaller sized matrix to test (like a 4x4). You could also submit an issue on their github.

Could this be a simple line ending issue? Are you using Windows to generate the file by any chance?

@terdon I am using Ubuntu, but good suggestion. Definitely been stung by the \r\n windows bug/feature before.

So I didn't find an answer,but I did find a way around if it anyone else has this problem. It does involve modifying the source code, but even if you don't understand c (I barely do) its still easy to do.
1) Download source code from github
https://github.com/cbcrg/tcoffee
Source code is here, view the dependency list - you will need all these libraries to compile the code for the change to take effect.
2) Add your custom matrix to source code.
Search matrices.h. In this version it is t_coffee/lib/matrices.h
You will see how they are defining their default matrices. Just add your matrix with the same format and give it a name. For example:
int customMT[]={
  4,
  0,  5,
 -3, -2, 17,
  0,  5, -3,  9,
  0,  0,  1,  1,  6,
 -2, -3, -3, -5, -4, 10,
  0,  0, -4, -1, -2, -3,  8,
 -2, -2, -5, -2,  0, -3, -3, 14,
  0, -2, -2, -4, -3,  0, -1, -2,  6,
  0,  0, -3,  0,  2, -1, -1, -2, -2,  4,
 -1, -1,  0, -1, -1,  2, -2, -1,  2, -2,  4,
  1, -2, -2, -3, -1, -2, -2,  2,  1,  2,  2,  6,
  0,  4, -1,  1, -1, -1,  0, -1,  0,  0, -2,  0,  8,
 -1, -2, -3, -1,  1, -4, -1,  1, -3,  1, -3, -4, -3, 11,
  1, -1, -2, -1,  2, -3, -2,  0, -2,  0, -2, -1, -1,  0,  8,
 -1, -2, -2, -1, -1, -1, -2, -1, -3,  1, -2,  0, -2, -1,  3,  8,
  1,  0, -2,  0,  0, -1,  0, -1, -1,  0, -2, -2,  0, -1, -1, -1,  4,
  1,  0, -2, -1, -2, -2, -2, -2,  0, -1,  0,  0,  1,  0,  0, -3,  2,  5,
  1, -2, -2, -2, -3,  1, -3, -3,  4, -2,  1,  0, -2, -4, -3, -1, -1,  1,  5,
 -5, -5, -2, -4, -1,  1,  1, -5, -3, -2, -2, -3, -7, -3, -1,  0, -3, -5, -3, 20,
  0, -1, -2, -1, -1, -1, -1, -1,  0,  0,  0,  0,  0, -1,  0, -1,  0,  0,  0, -2, -1,
 -4, -3, -6, -1, -2,  3, -3,  0, -1, -1,  3, -1, -4, -2, -1,  0, -2, -1,  1,  5, -1,  9,
  0,  0,  0,  0,  5, -4, -2,  0, -3,  1, -1, -1, -1,  0,  4,  0, -1, -1, -3, -1,  0, -2,  4};

Then, search for file tcoffee/lib/io_lib/io_func.c Find the function
int ** read_matrices( char *mat_name_in)

You will see a long list of if else clauses, add the name of your custom matrix using the same pattern as they do.
/*Read Header Matrices*/
if (strm(mat_name, "pam250mt"))vector=pam250mt;
else if (strm(mat_name, "idmat"))vector=idmat;
else if (strm(mat_name, "dna_idmat"))vector=idmat;

/*add custom matrix here*/
else if (strm(mat_name, "customMT"))vector=customMT;

3) Save the files, move into tcoffee/src/. Rebuild the files using the make file as described on the github installation. It is here you need make, g++, perl, and whatever else on the list as dependencies
make t_coffee

4) Now, align your file! The only thing here is to use the t_coffee of the source code you just edited, and not t_coffee that might already be on your path. So change t_coffee/src/t_coffee to whatever the relative path to this folder is.
tcoffee/src/t_coffee rhomboids.fasta -matrix customMT

And that's it!

Its a great question and presentation of the solution. It is an old last resort way of doing it, change the matrix and then re-compile. Reading the parsing code is the underlying issue. The only thing I would have suggested is get it working with the sample data first. It is a complicated table and if the parsing is sub-optimal and table not perfect it can trip the code. It only takes one part of the matrix to be missed and the whole thing falls over

Personally what I think has happened is the custom DNA matrix is probably heavily beta tested but the protein matrix much less so and I think thats what is being encountered here. Just to finally mention its also a good use of the custom distance matrix and interesting approach.


text:
How to remove sequence reads contains more than 2 X from multifasta file?
I have 5000 protein sequences in one multifasta file. I found more reads have gaps as X in their reads. So, want to eliminate those reads completely (Whole protein seq) from the file. I am keeping filter criteria as if a read contains morethan 2 X ( continuesly or anywhere in the read) should be removed. Thanks in advance for your help.
The input sequence looks like this
>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot2
ANSTVKKKKLLLYYYSSSEERXXXXXXXXXXXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

I want output Like this
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

Easy use Biopython "SequenceClean".You could do it as a Perl one liner or awk. Which of these solutions are you more familiar with?

Just for a record, X is not a gap but unknown amino acid.

Also, >Prot1  seems to have only one X, hence I suppose it should be kept as well, is that correct?

A Biopytnon solution that works for .gz files too :-). The script make a script filter_faa.py:
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') < 3:
            print(seq_record.format("fasta"))

Then usage is
python3 filter_faa.py my.faa > my_filtered.faa

The my_filtered.faa then will contain
>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.

it's generally better to explicitly close the open file handle, which won't happen for your gzip.open() case, could be modified, e.g. https://stackoverflow.com/a/52839332

Ha, I always wondered if there is an elegant soltion for opening it with/out gz. Cheers Chris!

It is an elegant solution even without closing the filehandle

Updated: as pointed out, this will not work for a multifasta as asked for in your original question. However it would work fine for a single line fasta

You could accomplish this with a single line of bash code if you would like (this assumes your fasta is in gzip format, if not just use cat):
zcat prot.fasta.gz | grep -E -v -B1 '(^>.*|^(.*X.*){3}$)' | grep -E -v '^--$'  - | gzip > prot.filtered.fasta.gzip

The -v option of grep inverts the regex, and the -B1 will print out the line occuring before the current context. The -B1 introduces -- between the matches, which is why the second grep is there. Then you can immediately pipe it back into a zipped file

I guess this would not work if the aa sequence had more than one line, or am I missing something?

You would be correct, In which case this would not work with OP's original question.  Grep would still be a useful tool for the filtering. You could probably filter out any sequences which havae More then 3 x's then filter out any Protein s that no longer have sequences


text:
Get a single number representing the contact probability between a pair of genomic loci using HiC data
Hello to the experts in analyzing chromosome structure data in HiC format,
I have a very basic question. I have a HiC file (specifics are mentioned below), using which I would like to obtain contact probability between a pair of genomic loci.
I obtained contact probabilities between the two loci at different resolutions available in the HiC file. However, rather than having a vector of contact probabilities, I would like to have a single number that appropriately represents the contact probability between the loci.
My guess is that there could be two ways to get a single number as the contact probability between two given genomic loci.

Should I take the highest contact probability?
or,
Depending on the lengths of the genomics loci of interest, may be the resolution that is close to that could be selected as a single number that represents the contact probability between the loci. Eg. with this strategy, if genomic loci of interest are ~5kb in length the contact probability at 5kb resolution is selected as a single number. Then, the issue is that the length of the genomic loci of interest could be lower than the lowest length of resolution available in the HiC file.
or,
Contact probabilities at different resolutions could be aggregated. Then I wonder what would be the appropriate way to aggregate contact probabilities, e.g. by multiplication or median or something else.

Just to reiterate, above these are my guess and I am not sure whats the best way.
So basically my question is: how can I get a single number that would represent the contact probability between a pair of genomic loci using HiC data?

Some specifics that might be useful in answering my question.

source of the HiC file: https://ftp.ncbi.nlm.nih.gov/geo/series/GSE63nnn/GSE63525/suppl/GSE63525_K562_combined_30.hic
organism: human (K562 cell line)
Resolutions in HiC data (bp): [2500000, 1000000, 500000, 250000, 100000, 50000, 25000, 10000, 5000]
Type of genomic loci: upstream non-coding regions of the genes
lengths of genomic loci of interest: average ~1000bp (range: 500-5000bp).
software I am using to "extract" data from HiC file: https://github.com/aidenlab/straw

The fundamental unit of a Hi-C contact is the nearest restriction enzyme cut site.
This is where each mate in a read pair from a Hi-C library comes from and is mapped to when converting sequencing data into a contact matrix.
Because the distance to the nearest cut site is variable at different genomic positions, and depth of coverage of any single locus is typically low, contacts are aggregated into bins and "normalized" to account for differences in cut site densities within those bins.
This normalization transforms the raw contact matrix into something resembling a stochastic transition matrix, where you can think of the elements of that matrix as "what is the probability that at any moment in time bin $i$ is in contact with bin $j$".
So the probability you're looking for is the $(i, j)$th element of the normalized contact matrix, where bins $i$ and $j$ contain your regions of interest.
If your regions straddle the bin boundaries, you may want to average those bins.
You get to pick the resolution of the contact matrix, but if you select a resolution that's too small, you run the risk of making your matrix too sparse and having no contact between those loci.
If you pick one that's too large, your matrix elements aren't really informative for your question at hand.
Given your 500-5000 bp range and median of 1 kbp, it may make sense to pick a resolution between 1 kbp and 5 kbp and use that.


text:
Identify differentially covered genes only between two samples
I have a question about finding differentially covered regions (coverage represents methylation level which goes from 0 to several thousands). I'm using enrichment based method which can be summarized with coverage per gene:
data <- matrix(sample(80), 20)
# Genes in rows
rownames(data) <- letters[1:20]
colnames(data) <- c("group_A_tr1", "group_A_tr2", "group_B_tr1", "group_B_tr2")

In data matrix each row represents a gene and each column represents a sample. There are two sample groups (A and B) with two technical replicates per each group. Problem is that we do not have any biological replicates.
My goal is to identify genes that are differentially methylated between two groups. I know that limma, edgeR, DESeq2 can be used in analysis like this, however I don't have enough samples. Basically I'll need to compare only two columns (after averaging technical replicates).
What method would be appropriate to work with data like this? Is it possible to treat technical replicates as biological ones?

Not really an answer but an extended comment... and most likely something you don't like to hear
I guess by technical replicates, it means taking the same biological sample and making 2 methylation libraries. If this is used as replicates in deseq2 or edger, the variance you are estimating is the technical variation that comes with preparing the library, bisulphite conversion etc..
Often the biological variance will be much larger. That is if you take 2 independent samples of group A, their variation in methylation is much larger than the technical variation you see from handling the samples. If the technical replicates are treated as biological replicates and differential analysis is performed, it will be seriously misleading.
I think both edgeR and DESeq2 doesn't estimate dispersion with 1 v 1 (which is right).
I would suggest maybe the following:

Run a fisher.test at every gene (methylated / unmethylated), get the odds ratio and p-value, rank this, but explain that the p-value is only for ranking, and will underestimate the variation.
A bit more refined, look for other studies that might be close to yours, get an estimate of the dispersion in their study and use a negative binomial distribution to again rank genes that show the biggest difference

You may run edgeR for methylation analysis without replicates (https://f1000research.com/articles/6-2055). But I also recommend you to have a look at the R DSS package (http://bioconductor.org/packages/release/bioc/html/DSS.html). It has a smoothing step which allows to modelize the biological variability by taking information from neighboring sites.


text:
Difference between PPI networks and Gene Co-expression networks
What exactly is the difference between protein-protein interaction (PPI) networks and gene co-expression networks? Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks, and so represent processed data. But could someone elaborate more on the differences?

Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks

Not really, no. PPI stands for protein-protein interaction. In PPI networks, each node is a protein and each edge (line) connecting it to another node represents an interaction. A very common example is that an edge represents a physical interaction between these two proteins. This means that there have been experiments to indicate that these two proteins come into contact in the cell and interact.
For example, they might form part of a protein complex, so their interaction was seen using chromatin immunoprecipitation, where the two proteins were captured together. Or, they may have been seen together using FRET, a technique that lets you detect when two molecules are in very close proximity.
Alternatively, some interactions can be theoretical, where in silico modelling has suggested that these two proteins can interact, even though we don't have any direct evidence that they do so.
Finally, yes, it is also possible to make a PPI network based on gene co-expression data, under the assumption that if two genes are co-expressed, then their protein products will be present in the same tissue. However, that doesn't mean they interact in any way, so it isn't really a very common way of building PPI networks, as far as I know.
In co-expression networks, on the other hand, each node is a gene and not its protein product. And here, the edges indicate the coexpression. Often the level of correlation between the two genes' expression data.
The most important difference between these two kinds of network, apart from the obvious mentioned above, is that they are used in different ways. A gene coexpression network lets you identify genes that are regulated in similar ways. It is a tool you can use to investigate the regulatory network of the cell, to see what genes are "turned on" or "off" together.
The PPI network, on the other hand, is a way of understanding the result of this gene co-regulation. You are not looking at how the gene network is regulated any more, but how the products of the genes interact with one another. So it is a very different focus.
You can maybe think of this in terms of school children and their parents. You could make a network where each node represents a pair of parents, and then draw lines between them if their children were born in the same year. That would give you a general overview of when this group of people procreated. However, that two children were born on the same day doesn't necessarily imply that the children themselves will ever come into contact. This could be the equivalent of the gene co-expression network.
Then, make another network where each node is a child and draw a line between any children who know each other. This would be the equivalent of the PPI network. You would likely find some correlations between the two, with children who know each other being more likely to have been born at around the same time, but the two networks give you very different information.

Co-expression networks are based on gene expression data. Genes are compared using their expression profiles over multiple samples. A common distance measurement is the Pearson Correlation Coefficient, assuming that co-expressed genes (genes expressed simultaneously in similar conditions) will have a high correlation coefficient.
In its simple form, a co-expression network is done like this. An expression matrix containing n genes x p samples is used in input (from microarrays or RNA-seq) to calculate a distance matrix of size n x n. Then the distance matrix is thresholded to keep only the most significant relationship (eg correlation >0.8). Gene pairs are connected if their correlation coefficient is above this threshold. Then all connections are visualized in a network, which can then be subjected to community detection algorithms to find more densely connected set of genes.
PPI are usually built on in vitro assays such as Yeast-2-Hybrid. But connections can also be inferred from literature data and expression data. You may have a look at the StringDB for example to see which kind of data can be used to connect proteins.

It is also possible to infer regulatory networks from gene expression data using specific data (eg mutants in transcription factors and controls) and specific metrics (eg bayesian inference). Infering a regulatory network is generally more challenging.


text:
How to split FASTQ reads without re-running `fastq-dump`?
I downloaded a few FASTQs from the SRA that were paired-end 76 bp.
When I look at the FASTQs, I get something like this:
@SRR7702334.1 1 length=152
AAGATNTTTTTCAGTATCAATCAGANGNNANANNCNNCTGCCCCAAATCCTGCAAAGAGCATGCAACGTCCTTAGAAGACAACGTCTTGCTCTGTCTCACCCAGGCTGGAGCATAGGGGCATGATCGATCATGACAAGGAGTTTACATTTTG

I forgot to add the --split-3 option to my fastq-dump command.
Is there an easy and fast way, without re-downloading the files, to split the reads in half and write them to two output FASTQs?
I tried the following using Biopython, and it'll work, but it's going to take a long time to finish.
from os.path import basename
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from gzip import open as gzopen
import argparse
from tqdm import tqdm

def split_read(read, pos=152 // 2):
    """
    Split read in half to produce 2 reads

    Parameters
    ----------
    read_id : Bio.SeqRecord.SeqRecord
        Read to be split
    """
    return {
        # R1/R2 as keys
        "R" + str(i + 1): SeqRecord(
            # first and second halves of sequence string
            seq=Seq(str(read.seq[(pos * i):(pos * (i + 1))]), IUPAC.ambiguous_dna),
            id=read.id,
            description="",
            letter_annotations={"phred_quality": read.letter_annotations["phred_quality"][(pos * i):(pos * (i + 1))]}
        ) for i in range(2)
    }

def split_fastq(fastq, prefix="split"):
    """
    Split reads from a FASTQ

    Parameters
    ----------
    fastq : str
        Path to input FASTQ file
    prefix : str
        Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)
    """
    # check for file format
    if fastq.endswith(".gz"):
        fq_handle = gzopen(fastq, "rt")
    else:
        fq_handle = open(fastq, "r")
    out_handles = {r: gzopen(prefix + "_" + r + ".fastq.gz", "wt") for r in ["R1", "R2"]}
    # load reads for random access
    records = SeqIO.parse(fq_handle, "fastq")
    for read in tqdm(records):
        # split read into each mate and fix header information
        mates = split_read(read)
        # write to each output file
        for m in mates:
            out_handles[m].write(mates[m].format("fastq"))
    fq_handle.close()
    for m in mates:
        out_handle[m].close()

if __name__ == "__main__":
    # parse command line arguments
    PARSER = argparse.ArgumentParser(description="Split a FASTQ into 2 FASTQs by splitting each read in half")
    PARSER.add_argument(
        "fq",
        type=str,
        help="Input FASTQ to split"
    )
    PARSER.add_argument(
        "-p",
        "--prefix",
        type=str,
        help="Prefix for output FASTQ files (`_R1.fastq.gz` and `_R2.fastq.gz` are appended)",
        default="split"
    )
    ARGS = PARSER.parse_args()
    # run splitting function
    split_fastq(ARGS.fq, ARGS.prefix)

```

oh... the sequences are joined.. this is nasty....i supposed you have a lot of files?

My comment above mentions the slow IO for the Python code (it was iterating at ~ 3600 reads/s, which, for a FASTQ with 500M reads would take ~1.6 d to complete), so I tried an alternative method that would cut down on the IO overhead. It's not super efficient, but I ended up using zcat | awk with a switch-case on the line number to stream the data.
# first pass to generate _R1 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 1, 76); break; case 3: print "+"; break; case 0: print substr($1, 1, 76)} }' | gzip -nc > {output}_R1.fastq.gz

# second pass to generate _R2 FASTQ
zcat {input} | awk '{switch (NR % 4) {case 1: print $1; break; case 2: print substr($1, 77, 76); break; case 3: print "+"; break; case 0: print substr($1, 77, 76)} }' | gzip -nc > {output}_R2.fastq.gz

I played around with trying to pipe each half of the combined read into its own file, but it was difficult to set that up and generalize it for the multiple FASTQs that I wanted to fix. I didn't feel like over-engineering something, and this was straightforward enough. It worked much quicker than I expected.
All in all this took ~4 h to complete, which means an iteration speed of ~ 70 000 reads/s, or about 20x faster than the Biopython solution.
If I engineered a solution that would have only passed through the FASTQ once, it would've taken half the time, but in either case it was likely shorter than re-downloading the entire FASTQ.
But I will probably make use of fasterq-dump next time, and make sure I specify to split the read mates.
I love how how extensive the default biopython fastq parser library is, but the slow IO for large sequencing files makes it unusable in a lot of cases, unfortunately. I wonder if there's a good way to use C bindings to make this faster, or whether you just need to use a different package/language.

I was thinking of two runs of awk, but I concluded that it can not be faster as you have to run it twice. Thank you for proving my intuition wrong.

I wish every bioinformatician knew this: never use the default biopython fastq parser on huge fastq files.


text:
Arlequin files not able to converge beyond 2000 steps for some .arp files
I am looking at 59 mitochondrial CO1 DNA sequences of length 563bp, but it is a haplotype file and it has 9 haplotypes within.
When I use .arp files in arlequin to run mismatch distribution it says that can’t converge beyond 2000 steps and also no graphs are produced. Can someone suggest what can I do to overcome this problem?

What you appear to have is a Bayesian MCMC for 2000 "steps" which isn't a lot and normally the term 'steps' is used for parsimony, so I'm a bit confused.
However, if a Bayesian analysis does not reach convergence, i.e. duplicate runs do not achieve the same results then the data space has not been solved so you get no result.
9 haplotypes does not sound much to me, it is possible there are too many sequences, it is also possible there is simply not enough variation (9 haplotypes ... how many alleles?).
Suggestions

Arlequin should provide a default convergence value which you can lower to force a result
Decrease the number of sequences but maintain the number of haplotypes
Try their example data set

Do you get any associated file such as posterior probability, any output file with a minus sign in front of it?


text:
Alternatives to phrap for merging contigs with overlapping ends?
I have bins from metagenomes generated with metabat2 on bitbucket/berkleylab/scr which are predicted as viral by MARVEL. MARVEL suggests to merge contigs with overlapping ends with PHRAP but the documentation on the tool doesn't seem great.
As it appears to be an assembler would it be okay to use another assembler? For example one of my bins has 4 contigs ~1000-5000bp in length. 3 of these seem to overlap with each other enitirely.

I found a tool MeGAMerge on github which might work. Has anyone used this?

you could try minimap2.  See here: https://github.com/lh3/minimap2/blob/master/cookbook.md#read-overlap

minimap2 seems like it will not work well for my contigs as they are not that long?

I found a tool called https://github.com/LANL-Bioinformatics/MeGAMerge that might work well? Anyone used it before?


text:
cellranger mkfastq with full path to --id flag
I have always used cellranger mkfastq to demultiplex 10x genomics runs manually, though recently the commands to do so have been incorporated into a script that should call the cellranger whenever it detects that a run is finished (and can be demultiplexed).
The main difference when calling manually compared to being called automatically by the script is from where the commands are being called.
So calling the commands manually from the directory where the samplesheet is and where the outputs shall be saved with:
/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=demux_10x --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc

does the job.
but now cellranger is being called (automatically) from /home/ngs. The error message that I am getting from cellranger is:
/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc
/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/bin
cellranger mkfastq (3.1.0)
Copyright (c) 2019 10x Genomics, Inc.  All rights reserved.
-------------------------------------------------------------------------------

/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/tenkit/bin/_mkfastq: line 531: __/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/.mro: No such file or directory
Martian Runtime - '3.1.0-v3.2.3'

Invalid name: /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ (only numbers, letters, dash, and underscore allowed)

According to the error message I cannot use slashes / in the --id flag.
Does anyone what I should pass to --id in order to get the outputs saved as
/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/?
I already know that if I just pass --id=demux_10x the outputs will be saved at /home/ngs/demux_10x. (again, I am calling the script that calls cellranger from /home/ngs/)
Obviously I could change to the desired output directory i.e. /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/ and do --id=demux_10, though I am trying to avoid that for now.

You should be able to use the --output-dir argument to specify a folder of your choice and use --id=demux_10x as before.

that works partially. All arguments passed to mkfastq were --id=demux_10x --output-dir=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc. The fastqs are being saved at /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10 while there is still a directory demux_10x and the file __demux_10x.mro at /home/ngs (from where the script was called). Hoping these files will be moved by completion of cellranger

The aforementioned files were not moved to the directory specified in --output-dir, even though the fastqs were saved there. On the _invocatoin file under /home/ngs/demux_10x there is an entry like interop_output_path   = null,, though apparently this cannot be specified by cellranger mkfastq. Guess the only solution will be to change directory back and forth

I remove quite a number of "such files" after my cellranger count runs, so this might be a deliberate design choice for some reason. BTW, the 10x customer service is very prompt, helped me in the course of a few hours in the past in case you would like to contact them.

I suspect the error message is a bit exaggerated.  Try removing the slash at the end of your --id

tried that way too though it did not work


text:
2x2 contingency for genomic feature enrichment
I have an enhancer list and I would like to see how it overlaps with histone lists from different cell types (different sizes) significantly. I would like to use a contingency table for testing but not sure which format:
Should I use 2xk format
                         Enhancer overlap       Enhancer not overlap
-Histone Cell Type 1             n                       n
-Histone Cell Type 2             n                       n
-Histone Cell Type 3             n                       n

or individual 2x2 formats
                                  Enhancer overlap      Enhancer not overlap
-Histone Cell Type 1 present           n                         n
-Histone Cell Type 1 absent            n                         n

                                 Enhancer overlap      Enhancer not overlap
-Histone Cell Type 2 present           n                         n
-Histone Cell Type 2 absent            n                         n

What do you exactly want to test? Maybe consider [edit] your question and clarifying the goal of this analysis.

the 2x2 is ok, you can test every histone (cell type) individually.... just control for multiple testing after

I think 2x2 is the easiest way. May be you are going to use fisher's exact test? Many tests that are designed for 2x2 cannot be applied to 2xk.

Some people would say that both schemes are inappropriate. The individual 2x2s are testing a slightly weird null, and the 2xk might be underpowered. If you're worried about this, I'd consider looking into multivariate categorical data analysis: https://journal.r-project.org/archive/2018/RJ-2018-015/RJ-2018-015.pdf. Contingency tables are still nice for looking at your data with human eyes, of course- for this purpose 2xk is more informative.


text:
Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) : cannot coerce class ‘"SeqFastaAA"’ to a data.frame
The example that I am trying to follow is this PGA tutorial.
I want to use the information in the analysis of the mzML raw files from the proteomics data analysis.
When I try to load the fasta file, I get the following error:
> msfile=openMSfile("/scratch/oknjav001/transcriptomics/proteogenomics/ms_rnaseqdata/msData/LTB/t16BCG.180731_C0_CH_T016BCG.mzML")
> msfile
Mass Spectrometry file handle.
Filename: 1807_C0_CH_T016BCG.mzML
Number of scans: 52332
> idfile <- runTandem(
+    spectra = msfile, fasta = dbfile, outdir = "./", cpu = 6,
+    enzyme = "[KR]|[X]", varmod = "15.994915@M", itol = 0.05,
+    fixmod = "57.021464@C", tol = 10, tolu = "ppm",
+    itolu = "Daltons", miss = 2, maxCharge = 8, ti = FALSE
+ )
Error in as.data.frame.default(x[[i]], optional = TRUE, stringsAsFactors = stringsAsFactors) :
  cannot coerce class '"SeqFastaAA"' to a data.frame

How can I go about fixing this error?
The example that I am trying to follow this PGA tutorial. I want to use this information in the analysis if the mzML raw files from the proteomics data analysis. When I try to load the fasta file I then get the above error.
msfile <- system.file("extdata/input", "pga.mgf",package="PGA")
idfile <- runTandem(spectra = msfile, fasta = dbfile, outdir = "./", cpu = 6,enzyme = "[KR]|[X]", varmod = "15.994915@M",itol = 0.05,fixmod = "57.021464@C", tol = 10, tolu = "ppm",itolu = "Daltons", miss = 2, maxCharge = 8, ti = FALSE)

sorry can you provide a bit more information, which package is this and what does your fasta file lloks like?

Hi Javan, welcome to Bioinformatics Stackexchange. To help you with your question, its often best to provide a minimal working example to see what causes your problem, and what you've attempted to try and solve it. You can edit your question to include code snippets, information about what libraries you're using, and what you expect your output to look like. At the moment, including the snapshot doesn't help us since we don't know what the input file is, what libraries you're using, etc. Not having that information makes it difficult for others to help answer your question

It looks like somewhere internally that spectra or fasta is expected to be a data.frame object, but is actually a SeqFastaAA object. The internal code is trying to convert is, but the class definition of SeqFastaAA has no method to be recast as a data.frame.
What is the dbfile you're using? Is is from the PGA tutorial you linked? And does your msfile coming from the mzML file need to be converted to a different object, like a data.frame before being used as an input for runTandem?


text:
How to convert BAM file to bigWig in Python?
I am manipulating some BAM files using pysam package which seems to be very fast and handy. However, I ran into problem when I am trying to generate bigWig files from BAM files in Python. I used bamCoverage from deeptools to do this previously, but it is a standalone tool that cannot be used from within my Python scripts.
Therefore, I am wondering what Python package would allow me to read in BAM files and convert them to bigWig files.
Thanks

deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:
from deeptools import writeBedGraph
from deeptools.getScaleFactor import get_scale_factor

wr = writeBedGraph.WriteBedGraph(options and input files)
func_args = {'scaleFactor': get_scale_factor(...options...)}
wr.run(writeBedGraph.scaleCoverage, func_args, output_file_name, ...)

Thanks, I am wondering what this scale factor means. Is it possible to calculate a RPKM or CPM like bamCoverage?

The scale factor dictates how the coverage track is transformed before writing the results to a bigWig file. Among other things, this then handles RPKM and other normalizations. Please notes that you'll have to pass an argparse object into it, which I understand is a bit annoying (I had't considered external use of that function).

pyranges:
import pyranges as pr
gr = pr.read_bam("your_file.bam")
gr.to_bigwig("out.bw", chromosome_sizes=pr.data.chromsizes()) # for hg19


text:
Getting weird plot with ggpubr package
I have a data frame for two groups and 9 variables
> head(aa)
  treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
> class(aa)
[1] "data.frame"
> str(aa)
'data.frame':   1026 obs. of  3 variables:
 $ treatment: Factor w/ 2 levels "post","pre": 2 2 2 2 2 2 2 2 2 2 ...
 $ variable : Factor w/ 9 levels "Sig1","Sig2",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ value    : num  0.117 0.189 0.219 0.128 0.171 ...
> tail(aa)
     treatment variable value
1021      post     Sig9     0
1022      post     Sig9     0
1023      post     Sig9     0
1024      post     Sig9     0
1025      post     Sig9     0
1026      post     Sig9     0
>

I am using this code
library(ggpubr)
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> my_comparisons <- c("pre","post")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y = 50)
Warning message:
Computation failed in `stat_signif()`:
missing value where TRUE/FALSE needed
>

But I am getting this weird plot

While something like this is my goal

I really don't know where I am doing wrong
I used ggplot where I got something like this but I need statistics that is why I have stuck on this

I then used this code
> DF <- aa %>% pivot_longer(.,cols = c(Sig1,Sig2,Sig3, Sig4,Sig5,Sig6,Sig7,Sig8,Sig9), names_to = "var", values_to = "val")
> View(DF)
> library(ggplot2)
> library(ggpubr)
> ggplot(DF, aes(x = treatment, y = val, fill = treatment))+
+     geom_boxplot()+
+     stat_compare_means(comparisons = c("pre","post")+
+     stat_compare_means()+
+     facet_grid(var~., scales = "free")
+

But I got another weird plot

But changing the scale I got this but gives one statistics

While I need one stat for each variable between pre and post groups
This solution gives to statistics
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment",facet.by = "variable")
> p + stat_compare_means(comparisons = my_comparisons)+ stat_compare_means(label.y =2)

When I edited my code
> p <- ggboxplot(aa, x = "variable", y = "value",color = "treatment")
> p + stat_compare_means(aes(group = variable))
>

But again one stat

The problem is the scale used: For the plot you called "weird" (first from the top), the scale is 50 and for the "ggplot only" (third from the top) the scale is 1.
You should play with the stat_compare_means(label.y = 50) bit, you can try setting the label.y parameter to 1.5 or 2.
EDIT BECAUSE OF THE ADDITIONAL REMARKS ON THE QUESTION:
Once again (and for sure the last time), I had to manually create example data as OP ignores repeated reminders on how to make a reproducible example with R.
library(data.table)
library(ggpubr)

text <- "
treatment variable     value
1       pre     Sig1 0.1173863
2       pre     Sig1 0.1888243
3       pre     Sig1 0.2191765
4       pre     Sig1 0.1277735
5       pre     Sig1 0.1711712
6       pre     Sig1 0.0000000
1       post     Sig1 0.1173863
2       post     Sig1 0.1888243
3       post     Sig1 0.2191765
4       post     Sig1 0.1277735
5       post     Sig1 0.1711712
6       post     Sig1 0.0000000
1021      pre     Sig9     0
1022      pre     Sig9     0
1023      pre     Sig9     0
1024      pre     Sig9     0
1025      pre     Sig9     0
1026      pre     Sig9     0
1021      post     Sig9     1
1022      post     Sig9     1
1023      post     Sig9     1
1024      post     Sig9     1
1025      post     Sig9     1
1026      post     Sig9     1
"
my_data <-fread(text)

> my_data
      V1 treatment variable     value
 1:    1       pre     Sig1 0.1173863
 2:    2       pre     Sig1 0.1888243
 3:    3       pre     Sig1 0.2191765
 4:    4       pre     Sig1 0.1277735
 5:    5       pre     Sig1 0.1711712
 6:    6       pre     Sig1 0.0000000
 7:    1      post     Sig1 0.1173863
 8:    2      post     Sig1 0.1888243
 9:    3      post     Sig1 0.2191765
10:    4      post     Sig1 0.1277735
11:    5      post     Sig1 0.1711712
12:    6      post     Sig1 0.0000000
13: 1021       pre     Sig9 0.0000000
14: 1022       pre     Sig9 0.0000000
15: 1023       pre     Sig9 0.0000000
16: 1024       pre     Sig9 0.0000000
17: 1025       pre     Sig9 0.0000000
18: 1026       pre     Sig9 0.0000000
19: 1021      post     Sig9 1.0000000
20: 1022      post     Sig9 1.0000000
21: 1023      post     Sig9 1.0000000
22: 1024      post     Sig9 1.0000000
23: 1025      post     Sig9 1.0000000
24: 1026      post     Sig9 1.0000000
      V1 treatment variable     value

ggboxplot(my_data,
          x = "treatment",
          y = "value",
          color = "treatment",
          facet.by = "variable") +
  stat_compare_means(label = "p.format")

And this code generates:

Thank you. now seems better but gives one statistics for all, what should I do to get one statistics in comparing each variable between two group?

You will need to play with the comparisons parameter within the stat_compare_means() function.

Sorry how? I have two groups pre and post which I placed them in comparison

A workaround would be doing facet.by = "variable" within your ggboxplot() call.

Thank you but gives no stat at all

Why do you call stat_compare_means() twice? And what happens when you use stat_compare_means(aes(group = variable))?

I removed extra stat_compare_means(), and stat_compare_means(aes(group = variable)) gives one stat

Thank you so much, are the stats are Kruskal-wallis?


text:
Using bioconductor from Python
Has anyone used bioconductor from Python? Is there any reason I would choose to use it from R instead of Python? It seems like there is a Python extension for it. Also, is there any reason to use bioconductor over Biopython?

don't quite get your question.. is there something you need to use one or the other for? The link you have uses rpy2, basically calling R from inside python.. not so easy for you to pass objects from rpy2 back into python

Bioconductor packages are R packages that can only be used directly from python via rpy2 or similar. If you've ever had reason to use rpy2 you'd probably realize quickly that while it works, it's no where near as easy as doing the same thing directly in R.
Biopython is a package for doing a number of useful things, but it can in no way come remotely close to performing even a percent of what Bioconductor packages do. That would be an impossible feat for a single python package.

Thanks for replying. So is there a way to work around not using R (and use Python instead) while still being able to do what Bioconductor can do? Is there an ecosystem of Python packages along with Biopython that can emulate what can be done with Bioconductor to a certain extent?

There is no python equivalent to Bioconductor. If you need to perform a statistical analysis you need to learn R.

Recently I have done something similar for my work, I am using sangerseqR which is s Bioconductor package in Python Environment using rpy2
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')

   #now using the sangerseqR package to the read the sequcence
   Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf")

I had some trouble installing rpy2 because of binary incompatibility issues. To solve this problem, I took the following steps:
First, I wrote the necessary functions in R. Then, I created a straightforward REST-API using library(plumber) and started the server on my local machine. On the Python side, I utilized requests to call the functions and receive responses.
I know that this solution is somewhat ugly, but it has works for me


text:
Evaluating amateur bioinformatics results
I've might have created a new algorithm for finding patterns in DNA. The technique has not been used before, and it managed to find a 33-mer with corresponding complement that exists both in E. Coli and in Salmonella. I think it's pretty neat but I don't know how to evaluate the results. In addition, I've developed this algorithm after a bioinformatics course where only the DNA sequences of E. Coli and Salmonella were provided, so I would like to test it on other sequences and have it evaluated by people who know what they are doing.
How would one get in contact with the right people for evaluating amateur bioinformatics result?

If you took a course, start by chatting with the Professor/Instructor who taught it. They will either be able to chat about things with you or put you in touch with someone locally.

It was one of those online Coursera courses. There are no real options to get in contact with the professor who created it.

What is the benefit of your approach? Is it faster than the myriad existing algorithms for pattern finding? Does it have better pattern-definition syntax? What kind of "patterns" are we talking about? What is this 33,er supposed to be? What makes it special?

I'm exploring the existance of superpermutations (of n = 4 = 33-mers) in DNA. I've created an array of all possible nucleotide superpermutations and I'm checking if sequences exist with a certain error margin. I've found a 33-mer that according to NCBI exists in 1839 unique organisms. I think this is interesting but I don't know if it means anything (with significance). The reason I'm exploring superpermutations is because they are the mother of all permutations (for which we look in DNA).

Wait, so you mean you just found some random sequence of 33nts whose only distinguishing feature is that you find it in many organisms? That doesn't seem strange, 33 is tiny, so you'll find loads of sequences that can be found all over the place at that length. Especially if you're not limiting to exact hits. Try generating some random 33nt sequences and searching for those, I bet many or even most of them will be found in hundreds of species.

Thank you. I was looking for an answer like that.

@terdon 33 is tiny, so you'll find loads of sequences I guess that depends how you look at it... with 33 nt you can make 4^33 ~ 7e19 combinations equivalent to about ~25 billion human genomes. If you pick one at random chances are it's never been seen anywhere (do I get it right?)

@dariober ah, but this isn't random at all. It's a superpermutation which I suspect will mean the sequence is relatively repetitive, making it far more likely to appear randomly in blast results. I admit I'm just speaking off the top of my head though (hence the comment instead of an answer). But yes, my assumption would be that finding a 33nt sequence present (as an imperfect match, mind you, as the OP described) in multiple genomes is not surprising per se. Finding exact matches in hundreds of genomes would be more intriguing, yes.

You need to understand the biological significance of the 33-mer locus.

I'd blast your 33-mer on NCBI and see what the biological significane is of the locus. If might be a transposon/ promoter region.
Look at the translation of the protein and if its part of a protein the location within the protein structure.
Blast the entire locus including 5' and 3' regions of the 33-mer target
Carefully assess the results for chance, e.g. conserved protein and the synonymous mutations just happened to drop out (Poison distribution).

Once you gain a greater understanding, it is likely to have been identified before and what you might contribute is a different method to obtain a known result.
In summary, you need the next level part of the course about translating a genetic correlation into a biological result/ mechanism. This isn't research grade work at present in my opinion.

Thank you for your answer, that NCBI is a really nice resource.


text:
Permanent error in installing R package
I want to install NNLM
install_github('linxihui/NNLM')

lang: error: unsupported option '-fopenmp'
make: *** [RcppExports.o] Error 1
ERROR: compilation failed for package ‘NNLM’
* removing ‘/Library/Frameworks/R.framework/Versions/3.6/Resources/library/NNLM’
Error: Failed to install 'NNLM' from GitHub:
  (converted from warning) installation of package ‘/var/folders/wk/qwgr9yqs6ts0czx__y_w8gch0000gq/T//RtmpA8PudK/filec5bce1f1309/NNLM_0.4.4.tar.gz’ had non-zero exit status
>
> install.packages('NNLM')
Warning in install.packages :
  package ‘NNLM’ is not available (for R version 3.6.3)
>

What should I do?
> session_info()
─ Session info ─────────────────────────────────────────────────────────────────
 setting  value
 version  R version 3.6.3 (2020-02-29)
 os       macOS Mojave 10.14.6
 system   x86_64, darwin15.6.0
 ui       RStudio
 language (EN)
 collate  en_GB.UTF-8
 ctype    en_GB.UTF-8
 tz       Europe/London
 date     2020-05-11

─ Packages ─────────────────────────────────────────────────────────────────────
 package      * version date       lib source
 abind          1.4-5   2016-07-21 [1] CRAN (R 3.6.0)
 assertthat     0.2.1   2019-03-21 [1] CRAN (R 3.6.0)
 backports      1.1.6   2020-04-05 [1] CRAN (R 3.6.3)
 broom          0.5.6   2020-04-20 [1] CRAN (R 3.6.3)
 callr          3.4.3   2020-03-28 [1] CRAN (R 3.6.2)
 car            3.0-7   2020-03-11 [1] CRAN (R 3.6.0)
 carData        3.0-3   2019-11-16 [1] CRAN (R 3.6.0)
 cellranger     1.1.0   2016-07-27 [1] CRAN (R 3.6.0)
 cli            2.0.2   2020-02-28 [1] CRAN (R 3.6.0)
 colorspace     1.4-1   2019-03-18 [1] CRAN (R 3.6.0)
 crayon         1.3.4   2017-09-16 [1] CRAN (R 3.6.0)
 curl           4.3     2019-12-02 [1] CRAN (R 3.6.0)
 data.table   * 1.12.8  2019-12-09 [1] CRAN (R 3.6.0)
 desc           1.2.0   2018-05-01 [1] CRAN (R 3.6.0)
 devtools     * 2.3.0   2020-04-10 [1] CRAN (R 3.6.3)
 digest         0.6.25  2020-02-23 [1] CRAN (R 3.6.0)
 dplyr          0.8.5   2020-03-07 [1] CRAN (R 3.6.0)
 ellipsis       0.3.0   2019-09-20 [1] CRAN (R 3.6.0)
 fansi          0.4.1   2020-01-08 [1] CRAN (R 3.6.0)
 farver         2.0.3   2020-01-16 [1] CRAN (R 3.6.0)
 forcats        0.5.0   2020-03-01 [1] CRAN (R 3.6.0)
 foreign        0.8-76  2020-03-03 [1] CRAN (R 3.6.0)
 fs             1.4.1   2020-04-04 [1] CRAN (R 3.6.3)
 generics       0.0.2   2018-11-29 [1] CRAN (R 3.6.0)
 ggplot2      * 3.3.0   2020-03-05 [1] CRAN (R 3.6.0)
 ggpubr       * 0.3.0   2020-05-04 [1] CRAN (R 3.6.3)
 ggsignif     * 0.6.0   2019-08-08 [1] CRAN (R 3.6.0)
 glue           1.4.0   2020-04-03 [1] CRAN (R 3.6.3)
 gtable         0.3.0   2019-03-25 [1] CRAN (R 3.6.0)
 haven          2.2.0   2019-11-08 [1] CRAN (R 3.6.0)
 hms            0.5.3   2020-01-08 [1] CRAN (R 3.6.0)
 labeling       0.3     2014-08-23 [1] CRAN (R 3.6.0)
 lattice        0.20-41 2020-04-02 [1] CRAN (R 3.6.3)
 lifecycle      0.2.0   2020-03-06 [1] CRAN (R 3.6.0)
 maftools     * 2.4.05  2020-05-10 [1] Github (PoisonAlien/maftools@6ddf526)
 magrittr       1.5     2014-11-22 [1] CRAN (R 3.6.0)
 Matrix         1.2-18  2019-11-27 [1] CRAN (R 3.6.3)
 memoise        1.1.0   2017-04-21 [1] CRAN (R 3.6.0)
 munsell        0.5.0   2018-06-12 [1] CRAN (R 3.6.0)
 nlme           3.1-147 2020-04-13 [1] CRAN (R 3.6.2)
 openxlsx       4.1.5   2020-05-06 [1] CRAN (R 3.6.2)
 pillar         1.4.4   2020-05-05 [1] CRAN (R 3.6.3)
 pkgbuild       1.0.8   2020-05-07 [1] CRAN (R 3.6.2)
 pkgconfig      2.0.3   2019-09-22 [1] CRAN (R 3.6.0)
 pkgload        1.0.2   2018-10-29 [1] CRAN (R 3.6.0)
 plyr           1.8.6   2020-03-03 [1] CRAN (R 3.6.0)
 prettyunits    1.1.1   2020-01-24 [1] CRAN (R 3.6.0)
 processx       3.4.2   2020-02-09 [1] CRAN (R 3.6.0)
 ps             1.3.3   2020-05-08 [1] CRAN (R 3.6.2)
 purrr          0.3.4   2020-04-17 [1] CRAN (R 3.6.2)
 R6             2.4.1   2019-11-12 [1] CRAN (R 3.6.0)
 RColorBrewer   1.1-2   2014-12-07 [1] CRAN (R 3.6.0)
 Rcpp           1.0.4.6 2020-04-09 [1] CRAN (R 3.6.3)
 readxl         1.3.1   2019-03-13 [1] CRAN (R 3.6.0)
 remotes      * 2.1.1   2020-02-15 [1] CRAN (R 3.6.0)
 reshape      * 0.8.8   2018-10-23 [1] CRAN (R 3.6.0)
 rio            0.5.16  2018-11-26 [1] CRAN (R 3.6.0)
 rlang          0.4.6   2020-05-02 [1] CRAN (R 3.6.2)
 rprojroot      1.3-2   2018-01-03 [1] CRAN (R 3.6.0)
 rstatix        0.5.0   2020-04-28 [1] CRAN (R 3.6.2)
 rstudioapi     0.11    2020-02-07 [1] CRAN (R 3.6.0)
 scales         1.1.0   2019-11-18 [1] CRAN (R 3.6.0)
 sessioninfo    1.1.1   2018-11-05 [1] CRAN (R 3.6.0)
 stringi        1.4.6   2020-02-17 [1] CRAN (R 3.6.0)
 survival       3.1-12  2020-04-10 [1] CRAN (R 3.6.3)
 testthat       2.3.2   2020-03-02 [1] CRAN (R 3.6.0)
 tibble         3.0.1   2020-04-20 [1] CRAN (R 3.6.3)
 tidyr          1.0.3   2020-05-07 [1] CRAN (R 3.6.2)
 tidyselect     1.0.0   2020-01-27 [1] CRAN (R 3.6.0)
 usethis      * 1.6.1   2020-04-29 [1] CRAN (R 3.6.2)
 vctrs          0.2.4   2020-03-10 [1] CRAN (R 3.6.0)
 withr          2.2.0   2020-04-20 [1] CRAN (R 3.6.2)
 zip            2.0.4   2019-09-01 [1] CRAN (R 3.6.0)

[1] /Library/Frameworks/R.framework/Versions/3.6/Resources/library

Update: I tried installing llvm, but don't know what put here to finish installing llvm-openmp pc-133-235:build fi1d18
$ cmake ../llvm -DLLVM_ENABLE_PROJECTS=openmp -DCMAKE_C_COMPILER=<C compiler> -DCMAKE_CXX_COMPILER=<C++ compiler>
-bash: syntax error near unexpected token `newline'

Answer from @atpoint, converted from comment:
There are dozens of threads in GitHub, SE and other support platforms that clearly explain what to do. (e.g. https://github.com/velocyto-team/velocyto.R/issues/2) The issue is that the compiler does not support openMP, so you have to install one that does.
Answer from @devon-ryan, converted from comment:
In conda on OSX we solve this by installing llvm-openmp, which provides the -fopenmp flag. Maybe the same is possible outside of conda (no clue). Don't compile it, find an easier way to install it.

try devtools::install_github('linxihui/NNLM'), if it doesn't work you need to install the gfortran binary on Mac.
If it still doesn't work, try to upgrade to 4.0.0.
When you have this kind of trouble, look at the issue page of the github page. You will almost always see that others have had the same issue.

Thank you installing gfortran did not help

Did you try to update R to 4.0.0 ?

This has nothing to do with binaries or the R version, it is a common issue with the gcc compiler on Mac which does not support openMP, see my comment above.


text:
umap and Louvain clustering on normalized data
I know that scaled data must be used for PCA for example as it is based on variance maximization.
However I'm wondering if it's the same case for UMAP ? If the data are single-cell RNA seq, after normalization can we do UMAP and Louvain clustering ?
I tried both and get similar results, however the Louvain clustering seems to be more adequate on normalized data than on scaled data.
By adequate I mean the clusters are the same but some are split into two, which makes sens looking at other results (like transcription factor analysis).

sorry what is the difference between scaled data and normalized data?

If the data are single-cell RNA seq, after normalization can we do
  UMAP and Louvain clustering ?

Typically people run PCA, UMAP and Louvain clustering on the normalised and log-transformed expression counts (but do marker gene and differential expression analysis on the non-normalised values). However, this remains controversial. I recommend reading "Current best practices in single‐cell RNA‐seq analysis: a tutorial"- it's a bit outdated, but still a good reference.

By adequate I mean the clusters are the same but some are split into
  two, which makes sens looking at other results (like transcription
  factor analysis).

The number of (or size of) clusters should not be used to interpret how good the clustering methodology is, since it is generally dependant on the input parameters, e.g. a higher resolution parameter leads to more clusters. Instead, display the clusters on the UMAP and see if they are continuous 'visually' and check what genes are differentially expressed between clusters and see if they make sense 'biologically'.

Since we both cited the same paper in our answers, do you think that a paper from 2019 can be called "outdated"? What new developments are not part of this overview?

The single cell field it moving so fast! e.g. Louvain clustering is now replaced by Leiden https://www.nature.com/articles/s41598-019-41695-z

I agree that there ara lot of developments happening fast (and as a matter of fact that Leiden works better than louvain). However, I think that the main points made by the paper, like the scaling question or how to normalize data are far more enduring and not addressed by any recent development.

Regarding the clustering evaluation bit above, there is plethora of metrics. I always keep an eye on the "average silhouette width", "adjusted rand index" (using SingleR annotations as labels) and also "clustering trees (easily obtained with clustree).

There is discussion if scaling the data (making the data range for all genes the same) is something you need for single cell data. The argument is that the expression differences between genes themselves are informative. Normalization on the other hand is always necessary.
The big single cell pipelines like Seurat or Monocle use both normalization and scaling as standard.
If you would like some of the clusters you get split into subclusters you could either increase the resolution parameter that the clustering functions usually have. Or you can take the cells from a cluster, then subcluster them and add the subclusters to the original clustering results.
For yor next question it would be great if you could provide the plots of the data because this way it is easier to understand what the problem is.

If you working on known cell culture, I recommend using classification approach instead of clustering. One library that I used to publish my paper is SingleR. I checked the expression x cell type matrix by both approaches (Seurat Clustering and SingleR), the classification approach make much more sense.


text:
Extracting specific sections from FASTQ determined by positions
Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions.
For example, say I want to extract the nucleotides at position 1, 3, 5-7 from each sequence:
Input:
@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

Output:
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)
Ideally, a command should work this way:
subselect 1,3,5-7 input.fq > output.fq

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5` and pasting) but maybe out of an box solution does already exist which evaded me.

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:
$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.
IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):
#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:
perl foo.pl file.fq 1,3,5-7

This looks like a job for samtools fqidx (http://www.htslib.org/doc/samtools-fqidx.html), assuming that your file is not enormous... read the warning in the link above.
usage would be
samtools fqidx yourfile.fq <region>
where region is specified using the fastq record identifier and coordinates.

How would a region like what the OP wants be written? You can run samtools fqidx file.fq id1 to get the data for id1, but how would you specify only the 1st, 3rd, 5th, 6th and 7th residue as requested by the OP? The docs don't show the format. And wouldn't this need each id in the fastq file to be mentioned individually? How would you use this for what the OP asked for?

I second terdons, fqidx is certainly a useful tool, but it doesn't address the specific question here.


text:
Should biological replicates be the most similar pairs in an RNAseq experiment?
The attached figure (from deseq2) shows the sample to sample distances of a RNAseq experiment with 4 conditions (A,B,C,D) at 2 time points (0h,4h) with 2 biological replicates.
I am a little bit puzzled that the biological replicates are not always the most similar pairs. In a perfect world biological replicates should be identical but I assumed that at least they should be very similar.
Is this something I have to worry about, or is this more or less normal?

This depends on the magnitude of biological variation and technical variation in your data.

In summary you basically have noisy data and to simply answer your question it depends on the genetic diversity of the samples you are dealing with. Note, from my recollection this was confusing to start with because the "replicates" column was labeled as "treatment", but I could be wrong.

Revised answer, firstly there is zero stationarity - thus its a straight time-series.
Summary As a time series

The statistical noise between A and D is lower than the difference in phenotypic expression between each and all othe replicates in the analysis ABCD - FOR DAY 0 VERSUS DAY
4, so they cluster and are distinct.
The statistical noise between B and C is greater than the phenotypic variation between their expression profiles for DAY 0 to DAY 4. I suspect they are genetically very similar. Thats the conclusion, you can spice it up with different analysis, but I'd leave it because might find a lack of statistical significance in your core observations (A and D).
Within the time-points the data is very noisy and doesn't appear to show a distinct signal. Thats okay, this is essentially a time series.
Sample D is distinct in all scenarios (which is good)

Result To reiterate and expand a bit, the result you have is that sample A and D biological replicates are consisent within the analysis you are performing, i.e. day 0 verses day 4 (below). Thus they cluster together and this is proof of distinct phenotypes above the statistical noise (A with A and D with D). An intra-time point dendrogram would have a different result. The replicates samples B and C are not consistent and the interpretation is that the phenotypic expression profiles are so similar they are essentially identical within the boundaries of statistical noise.
Analysis The dendrogram is combines day 0 to day 4 analysis with day 0 to day 0 analysis and day 4 to day 4 analysis. HOWEVER, the clustering algorithm will correct for this and separate based on the largest distance (day 0 versus day 4). The signal that you get in your dendrogram is the largest signal in the data set which is day 0 vs day 4, and this is how UPGMA clustering works (I probably have previous posts on the statistical algorithm of UPGMA). So actually that is fine.
Advice With data this noisy I'd go no further in attempting to deploy critical probablity (a p value)  because this data could be the result of a complex distribution. This would require a statistician. Personally
 I'd upload the PCA as a separate post, but its just my opinion, because UPGMA versus eigenvalue based PCA are statistically quite different approaches.

Thank you. X-axis is just the same as the y (flipped). This is a all-versus-all pairwise comparison.

Hang on, none of the conditions are replications. Each experiment is different?

It doesn't look bad, the clustering does not seem to be driven by batch (replicate). time has a strong effect and it is also clear that D has rather distinct biology than others. Do you also have a PCA or an MDS plot? Such a plot where condition is represented as colors and time as shapes for example can help getting the full picture? And if your design matrix enables, you can add all of these variables in your formula to control for confounders.

I think we need to know biologically what 4hours representsds, confirm this is not a time series and we could could expect 'stationarity'

Yes I have a PCA as well, I will upload later today. PCA shows the same effect, very strong separation by time.
But I am not so interested in this time effect, rather I would compare A vs B for 0h or for 4h.
For deseq I used the simple all versus all option (galaxy), so galaxy does is only aware of a condition that is the concatenation of condition+time.

Okay I get it, I've updated.