text:
Joint variant calling of two sample sets with different coverages (~5x vs ~25x)
Let's say 50 samples were sequenced at 25x. After a while, the group realizes that the experiment is under-powered and adds 100 samples more, but at a coverage of 5x. I'm talking about WGS.
The scientist runs the samples following GATK's best practices until he/she is done running haplotypeCaller on them.
Then the question is how to proceed:
a) combine all 150 gVCFs and do joint calling. Then do site filtering followed by genotype filtering.
b) Do joint calling with the high coverage samples and in parallele do joint calling with the low coverage samples. Then do site filtering, merge both VCFs and filter by genotype.
c) combine all 150 gVCFs and do joint calling. Split VCF into two according to coverage and do site filtering. Merge both VCFs and filter by genotype.
Option "a" sticks to GATK's recommendations, but it ignores the high difference in coverage between sample sets. The opposite is the case for option "b". And option "c" is like a hybrid approach of "a" and "b".
I guess what I want to ask if it is an issue to do joint-calling and then site-filtering when groups of samples have large differences in coverage.

I’m guessing a bit here, but I don’t think differences in coverage should matter, since coverage is accounted for during the joint calling stage (by looking at the read support for each genotype). So I’d guess b would be a better option. Of course, there may be downstream effects of having one batch at higher coverage than the other, but you should be able to account for that in the filtering stage.


text:
workflow for processed RNAseq data analysis in R
I'm learning to use R in data analysis. I'm getting fluent in baseR and the tidyverse, but thus far only dealt with medium throughput plate-based experiments.
I am currently trying to learn how to work with a data set of normalized RNAseq results. I have a matrix consisting of tens of thousands of rows corresponding to genes. The first two columns identify 'gene_ID' and 'gene_name'. The remaining 96 columns are different samples (different treatments, timepoints, replicates, ...).
I would like to analyze this dataset to identify differentially expressed gene based on the different conditions, perform some cluster analysis, run GESA on these clusters, and generate heatmaps at various stages of the process.
Can anybody recommend the best libraries to look into in order to achieve these goals? I started looking at several tutorials, but much of them deal with the earlier steps of RNAseq data analysis and don't really focus on these later steps once the data is already processed and cleaned up.

What tutorials have you read? DESeq2 starts from raw RNAseq results and goes into DE analysis, after which you can look at GSEA and similar approaches. As for data types/containers for your data based on its size, you can use plain data frames or SummarizedExperiment objects.

I have not yet worked myself through any specific tutorial since there seemed to be so many tools available. I started off with my 'standard approach' of turning everything into long data frames with columns for the various treatment types, but that became completely unworkable (8 columns with 5.5 million rows). So I was hoping to get a pointer about which package(s) and data formats I should start to investigate first.

If y have normalized data, I think limma is your best bet for finding DE genes.  EdgeR and DESeq are the most commonly used, but they require raw reads counts.

Thanks @swbarnes2. I'll start my investigations there.

source: Stack Exchange

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

sort: votes


text:
How do I force Cromwell / WDL to put input files in the same directory?
I have a cromwell task that, among other variables, takes these three variables:
   File Ref             # Reference genome fasta
   File RefFai          # Reference Genome index
   File RefDict         # Reference Dictionary

While all three files make it into the task's "inputs" directory, they end up in different sub-directories of inputs:
ME $ find inputs
inputs

inputs/-181926642
inputs/-181926642/Sample.vcf.bgz.tbi
inputs/-181926642/Sample.vcf.bgz
inputs/-1548908437
inputs/-1548908437/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.dict
inputs/44588251
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa
inputs/44588251/GCA_000001405.15_GRCh38_no_alt_short_headers_nonACTG_to_N.fa.fai

Which means that when I call GATK and pass it the path to the .fa file, it gets rather cranky that the .dict file isn't in the same directory.
How do I force cromwell to put all three files in the same directory?

We did have a similar problem and there is a solution to it... however, I have hard time recalling it, sinc eit was a while ago (and we abandoned WDL). Perhaps you could paste the code for the task, which would make it easier to see, if anything is missing.

@Vadim the code does not reference RefFai or RefDict, only Ref.

I am talking about the wdl code: one must pass all the files explicitly to the task, even if they do not appear in the command line. Every task is run from a separate folder.

My hacktastic fix for this is to use basename to get the name of each file, then do
refName=`basename "${Ref}"`
ln -s "${Ref}" "$refName"

for each file.  This puts aliases to each file in the executions dir, and then I pass "refName" to GATK rather than "${Ref}".
Not a good solution, but a solution

as for wdl version 1.0 (other versions untested)
you can pack your reference and related files into a struct:
struct Reference {
    File fasta
    File dict
    File fai
}

You will pass the Reference to your tasks instead of File.
Inside your workflow/config define:
Ref example = {
    "fasta" : "reference.fasta",
    "fai" : "reference.fasta.fai",
    "dict" : "reference.dict"
}

Inside the task you just use the file want use by example.fasta.

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


text:
How to check if a given gene is expressed in a group of microarray samples if I do not have control group to compare with?
I have a microarray gene expression dataset consisting of placenta samples. I want to check whether these placenta samples are mixed with maternal decidua tissue. I have marker genes for decidua samples SCARA5, PIGF, SLPI. How to check whether these genes are expressed in placenta samples? In other words, how to check if a gene is up or down regulated without control group?

Which microarray? Some have negative controls that can be used to assess the probability that a gene/transcript is expressed.

I look at theese datasets: GSE73374, GSE73685, GSE9984, GSE122214, GSE22490, GSE35574, GSE37901. The microarray platforms are Affymetix and Illumina: affy_hugene_2_0_st_v1, affy_hugene_1_0_st_v1, affy_hg_u133_plus_2, illumina_humanwg_6_v2. Unfortunately I do not have access to illumina raw data, only to the preprocessed one. Can I still get negative controls from it?

At least some of those have the .CEL files included, so you have all of the data.

@DevonRyan I've looked .CELL files' description and indeed there are control probes. Could you give me a hint, how do I move from these raw data to the lists of up- and down-regulated genes?

Finding up and down regulated genes requires using R packages like affy and limma. Please search for tutorials.

Thank you @DevonRyan I will look them up more closely. But aren't these methods require having case and control groups of samples and not just one group like in my case?

Your previous comment asked for up and down-regulated genes. That requires some sort of treatment groups, yes.

Oh, I see now. I said up- and down-regulated, but I meant expressed and unexpressed relative to control probes. My bad.

You ask about which "genes are expressed" and then you mention "if a gene is up or down regulated". These are different, and given your application I think what you actually want to know if these marker genes are expressed. This is not a question of up/down-regulation relative to a control, so you do not need control samples for this.
It seems most of your datasets are Affymetrix, and there are R/Bioconductor packages that will give you the probability a gene is expressed. One such I found is frma (https://bioconductor.org/packages/release/bioc/html/frma.html), which returns an object "weight" that is "a vector of weights which roughly correspond to the probability of expression for each gene." I haven't worked with microarrays in a while, but I recall many other R/Bioconductor packages offering this functionality.
A detection p-value is calculated for Illumina array probes in Illumina's native BeadStudio software.


text:
How to call this : some dots in volcano plot aligned a line after add a very small value to padj
when I plotting the volcano plot, I got some dots very top of the Y-axis because their padj nearly reached 0（I checked the data, I think they are reliable）, so I followed the advice to add a very small value to padj(machine$double.xmin), but still got those dots formed a line. I think I heard these lines have a name like Texas arm(?) or Tech arm(?) something before, but I can't remember. Can someone tell me what's the name?
Please help me. here is the plot :

If there's a name for that then it's so obscure that if you used it in a paper the reviewers would have to google for it.

@DevonRyan Yes, I think you are right. Maybe I'm misremembering. Anyway, thank you for your help!


text:
Multi Omics comparison RNA seq : Condition specific gene signature
The authors of an article are comparing different cell type and the organ type and they show the ontology enrichment of the RNA-seq:

Cell-type-specific and organ-specific marker genes were identified
using a two-step procedure. First, we performed pairwise differential
expression analysis for each of the three type of structural cells,
comparing them across organs. Second, based on the resulting pairwise
comparisons, we identified for each organ those genes that were
upregulated in each cell type compared to at least five other organs.
Pairwise comparisons between organs were performed using the limma
package in R, separately for each of the three structural cell types.
Significantly differential genes were selected based on statistical
significance (adjusted P < 0.05), average expression (log2(CPM) > −1)
and sequencing coverage (median number of reads greater than 10 in the
group with stronger signal). On the basis of these pairwise
comparisons, we counted the total number of times each gene was
upregulated in a specific organ compared to all other organs. Genes
that were upregulated in comparison to five or more other organs were
selected as marker genes of the corresponding organ.

Now how is this part done ? "Genes that were upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ.*"
The figure for this results is this:

What part of "Genes that were upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ." is unclear?

yes sort of .this part "upregulated in comparison to five or more other organs were selected as marker genes of the corresponding organ" , This one  "we identified for each organ those genes that were upregulated in each cell type compared to at least five other organs." so if i do pairwise comparison then i would be getting upregulated in one vs the other. So here how do they show enrichment to each individual cell type such as brain liver etc etc and from the same data they go for enrichment ..

im curious how they do to get single cell type gene expression which is Upregulated ..for example if i have 3 cell type either will do pairwise or will do across the grroup such as LRT in deseq2. if i do pairwsie i will always get UP or DOWN with respect to my control isnt it?

From their methods section:

On the basis of these lists of marker genes and receptor–ligand pairs, we inferred potential cell–cell interactions for all pairs of one structural cell type and one haematopoietic immune cell type, quantifying the enrichment for known receptor–ligand pairs among all pairs of marker genes between the structural cell type and the haematopoietic cell type (Fig. 2a). First, we counted all pairs of marker genes for each pair of structural and immune cell types. Second, we calculated the fraction of these gene pairs that were annotated as receptor–ligand pairs. Third, we tested whether this fraction was greater than the fraction of annotated receptor–ligand pairs across all pairs of genes. Fisher’s exact test was used to obtain P values and odds ratios, as implemented by the function ‘fisher.test’ in R. P values were adjusted for multiple testing using the Benjamini–Hochberg method, as implemented the function ‘p.adjust’ in R. Finally, significantly enriched pairs (adjusted P < 0.05) of structural and immune cell types were connected by edges to generate a graph of cell–cell interactions.

In addition to testing for significant enrichment (as described in the previous paragraph), we calculated aggregated receptor and ligand gene signatures (Fig. 2c) based on a manually curated list of immune-related receptors and ligands (Supplementary Table 4). These gene signatures were derived by normalizing expression values for each gene and averaging across all genes in a given set of biologically related receptors and ligands.

The actual values plotted then seem to be from EnrichR, since they mention that elsewhere in the methods.

"The actual values plotted then seem to be from EnrichR, since they mention that elsewhere in the methods." Okay...so they simply taking this "These gene signatures were derived by normalizing expression values for each gene and averaging across all genes in a given set of biologically related receptors and ligands. , if that is the case then the the bubble plot shouldn;t be similar?


text:
FASTQ file reads extraction paired ends
I have paired end reads and want to extract only those which have this sequence within them:
TGTATGTAAACTTCCGACTTCAACTGTA

I tried using grep -B1 -A2 "TGTATGTAAACTTCCGACTTCAACTGTA" file_1.fq | grep -v "^--$"
and the same from the mate pair. After this I tried aligning them , but they dont align because the reads headers are different.
Can someone help in extracting the mate pair of one read from the other file (2.fq), so that they can be aligned ?
regards.

Can you sit down and write out exactly what you want to do?  I think it's pretty clear why what you did is not returning matched pairs of reads, but it's not clear to me if you see why.

As I wrote, that I have paired end reads. I want to extract reads from both pairs which have in inserted sequence TGTATGTAAACTTCCGACTTCAACTGTA in them. It belongs to a transposon along with guideDNA. The fastq files are stored as forward and reverse reads. When I executed the script it brings out the reads which have insertion within them. However, they dont align with Bowtie2 as the reads are not from the same pairs.

What I'm trying to get at is, do you want reads which lack the sequence if their mate has it?  Because if the answer is "no", why are you trying to align paired reads when you are throwing away a member of a pair?  And if the answer is "yes", you are clearly throwing those reads away.

Your question can be rewritten as, "how do I get the union of reads in two files?" Since you have already subset the mates, the steps are:

Process each file, storing the list of reads names in each.
Take the union of these.
Iterate through the original/non-subset files and write only the reads with names in the union to the output.

The files will be both subset correctly and synced properly such that bowtie2 or any other aligner can correctly process them.

Hi Devon....is there a tool which can do this ? Kindly let me know.

Probably, but in the time it would take to search for it it'd be faster to write the 10-20 lines of python (or whatever you prefer) to just do it.

If you want to stick with shell scripting: comm -12 file1ids file2ds should give the union, followed by a similar grep like you used but with -f union.txt to give the patterns in a file. The undocumented --no-group-separator argument to grep may help you too if you're using GNU grep.

Alternately, you could go through both files simultaneously, and output both reads if either of them matches.

Based on my reading of the question, I'd guess they're interested in the intersection of reads, rather than the union

@gringer Could be, I had used the intersection originally in my reply but then figured OP wanted reads with the pattern in either mate. We'd need clarification I guess.

It's a little bit unclear what you're asking, but if you want to repair paired-end files that have become disordered (or the linking has been lost), then the repair function of bbtools should work. This will pair up any coincident read IDs, and create a separate file for any reads that have no pairs in the other file:
https://jgi.doe.gov/data-and-tools/bbtools/bb-tools-user-guide/repair-guide/
repair.sh in1=broken1.fq in2=broken2 out1=fixed1.fq out2=fixed2.fq outs=singletons.fq repair

If grep is used on one file, and repair used on that file in combination with the original file from the matching read, then the relevant read pairs should appear in the output:
repair.sh in1=filtered_1.fq in2=file_2 out1=fixed_1.fq out2=fixed_2.fq outs=unpaired_2.fq repair

It's not clear to me that both halves of any pairs will be caught by that grep, so there won't be anything repairable.

That's the "unclear what you're asking" bit.


text:
How to extract reads with INDELs > a given size?
I'm trying to modify this https://www.biostars.org/p/253774/
To get reads with deletions > 20bp
I think this gives reads with exactly 20bp dels:
samtools view -h in.bam chr1:12340-200000 \ | awk '$1 ~ "^@" || $6 ! "20D"' \ | samtools view -b - > out.bam

is it possibe to do >?
Thanks!

There's no simple way to do that with awk, since you would need to parse the CIGAR string and iterate over the operations. In python using pysam you would iterate over the reads and filter as follows:
keep = False
for OP, LEN in b.cigartuples:
    if OP == 2 and LEN > 20:
        keep = True
        break

if keep:
    ...do stuff...


text:
What is "unmapped read segments" in the output of samtools idxstats?
samtools idxstats produces a four column output  (see here)
The output is TAB-delimited with each line consisting of reference sequence name, sequence length,
# mapped read-segments and # unmapped read-segments. It is written to stdout. Note this may count reads
multiple times if they are mapped more than once or in multiple fragments.

What is the meaning of unmapped reads here? Naively, I would take it to mean all the reads that didn't map to a specific feature... however in a sample with millions of reads most of them should be counted as such, which is not the case in my calculations, where the number of unmapped reads is of the same order as the number of mapped reads.
I will very much appreciate explanations from people experience with samtools.

An unmapped read segment is any entry with bit 4 of the flag field set (i.e., it's unmapped). Assuming the reference to which you aligned reasonably matches the sample from which the reads were generated, less than ~10% of to total reads should be unmapped.
Please note that BAM files have no idea what a "feature" is. They only contain information about reads and the chromosomes/contigs to which they align.

I have mapped to a gene catalogue, i.e., I have many references, and every line in the idxstats output correspondes to one such reference. I understand how a read mapped to a reference/gene/feature is associated with it, but a read unmapped is unmapped to all of them...
But what you wrote about flags makes a lot of sense - I have to look what happened at the mapping step. Thanks.

For paired-end reads, some aligners will assign an unmapped read to the chromosome of its mate if its mate aligned (the unmapped bit in the flag will still be set).


text:
Getting System.IO.FileNotFoundException with MaxQuant running on Linux with mono
I am getting this error with MaxQuant 1.6.14 running with mono (5.14.0.177=he1b5a44_1001).

System.IO.FileNotFoundException: Could not find file "/home/ubuntu/run-dir/combined/combinedRunInfo"

And I not sure why I am getting this error. The process has read and write permissions to that location.

...seems the mqpar file was not compatible with the maxquant version.
Make sure the mqpar.xml file shows the same version number as the MaxQuant binary.
MaxQuant is not really made for the command line. So, the errors shown are sometimes quite misleading.


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

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

Permutation: free
Number of permutations: 5000

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

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

Many thanks

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

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

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

This is done; thanks!

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

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

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

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


text:
How to rename Read Groups
I periodically run into problems with Read Groups. I am in a desperate need for a tool that would allow me to fix them.
This time, when I was mapping reads, I formatted all the information in the Read Groups right, but the name of the read groups (which I, unfortunately, chose to be RUN_LANE, not realizing that some of the samples were on the same lane in the same run). And once I tried to do a joint genotyping, freebayes got mad at me.
So, is there an easy way to rename read groups?
Piccard tools' AddOrReplaceReadGroups won't do the job as some of the samples were multiplexed on different lanes/runs and this information should be preserved.
Basically, I need all bam files with read groups like this
@RG     ID:HJ5FNDSXY_L2 SM:WW5-6        PL:ILLUMINA     LB:WW5-6
@RG     ID:HJ5JVDSXY_L2 SM:WW5-6        PL:ILLUMINA     LB:WW5-6
@RG     ID:HJHGLDSXY_L4 SM:WW5-6        PL:ILLUMINA     LB:WW5-6

to change to something like (so the ID string will be unique for every combination of sample, run and library)
@RG     ID:WW5-6_HJ5FNDSXY_L2 SM:WW5-6        PL:ILLUMINA     LB:WW5-6
@RG     ID:WW5-6_HJ5JVDSXY_L2 SM:WW5-6        PL:ILLUMINA     LB:WW5-6
@RG     ID:WW5-6_HJHGLDSXY_L4 SM:WW5-6        PL:ILLUMINA     LB:WW5-6

There is a similar question asking about adjustments of read groups using pysam. That is actually quite close to what I need, but it's rather hacky solution I would rather avoid if I can.
I also vaguely remember seeing a tool for fixing read groups problems while ago on Twitter, but I can not find it now and I don't know if I just made it up or not.

The most annoying part of this problem is that the lowest effort solution is to remap the reads with corrected RG IDs (but I always feel bad for burning computational hours in vein).


text:
how to repeat records in fastq n times efficiently?
How to iterate/repeat a record n times in a fastq file using bioawk? I wrote a python code using biopython, but it is very very slow. So, I am wondering if I can get some help by using bioawk. Thank you so much!
for example, in the file.fq file, I'd like to repeat read_id1 for 100times, read_id2 for 100times, ..., and write them to a new file.

Hi and welcom to BS, could you please a little bit more what you problem you want to solve? Do you want only specific record or all records of the file? Maybe you can show your python code for clarity.

You might want to show your biopython implementation. While biopython can be slow, your solution may be sub-optimal

There are quite a few people complaining about the speed of Biopython... I think pysam is much faster.

Not sure why you would want to do this, but if the aim is to repeat all records in the fastq/fasta files n times, the most efficient way I can think of is just to concatenate the orginal file n times, similar to this thread.
i.e adapting the accepted answer would look something like this (I have not tested this):
for i in {1..n};do cat file.fq >> file_ntimes.fq ; done

Update based on comments:
To replicate each fastq entry back to back, maybe this hacky unix command can be adjusted (say for n=4):
n=10
paste - - - - <file.fq | awk '{for(i=0;i<$n;i++)print}' | tr '\t' '\n' > file_10times.fq

The idea is to paste blocks of 4 lines (i.e. 1 fastq entry) into one line and replicate this line n times. tr is reversing the paste action (since the separation by paste is a tab character). This won't work if fastq file for some reason has tabs in it which I don't expect but didn't test either.
Note, for fasta file this would be paste - - rather than paste - - - -
I haven't tested this, but if it works for you I could polish it up if needed.
Update2
Based on Xioafeis comment this is solution seems to work as well:
bioawk -cfastx '{for(i=1; i<=10; i++) {print "@"$name"\n"$seq"\n+\n"$qual}}' file.fq > file_10times.fq

Thank you so much for your reply!
But I don't think I made my question clear. I want to repeat rd_1 for 100 times, rd_2 for 100 times, ... in this order, not (rd_1, rd_2, ...) x 100 times.

My goal is to pairwise each of read record in file_1.fq with each of record in file_2.fq into 2 new fastq files and get overlap between them. Actually, for the one that you suggested that is what I want to for my file_new_2.fq, which I have already figured out. But not the one that I am asking. I will edit my question with a figure. Really appreciate for your reply!

I did wrote a python code using biopython, but it is very slow I guess it slurped whole file into dictionary, so very very slow and seems running forever

Thank you so much! I figured it out by using bioawk.

bioawk -cfastx '{for(i=1; i<=10; i++) {print "@"$name"\n"$seq"\n+\n"$qual}}' test_repeat.fq > test_repeat_10times.fq

I see, did my suggested solution work as well? Feel free to upvote or accept the answer in that case! Also Happy to add you solution as 3rd option if you want to have it as offical answer?


text:
How to get matching pattern along with ID in a single command in grep?
I have a file containing multiple fast sequneces. For a specific consensus pattern as input, I extracted all the matching patterns from target fasta sequences with
grep -o -E "CC[GT]AAA[GC][AC]TT[GC]" input.fasta

However, the above command will retrieve just the matching sequences from fasta sequences and I  also wanted to get the corresponding fasta header of each match.
For example, if input.fasta  file is something like this,
>Gene 1
TGATGAAAAATGATAGAT
ATTGGGGGAAAAAAAAAT

>Gene 2
TTTCCTAAAGATTGT
AAATTTAAAAATGTTTTT

(Gene 2 has matching subsequence CCTAAAGATTG)
Output:
CCTAAAGATTG   Gene2

I prefer a solution with grep. But other possible solutions also helpful.

Use this Perl one-liner:
perl -lne '$id = $1 if /^>(.+)/; ($m) = /(CC[GT]AAA[GC][AC]TT[GC])/; print join "\t", $id, $m if $m;' input.fasta

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-n : Loop over the input one line at a time, assigning it to $_ by default.
-l : Strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches


text:
Peak Specific to condition in ATAC seq data
How it is define in terms if accessibility. For example in expression data when I compare A vs B either A or B as reference we define that genes are getting differentially expressed to A or B with our choice of cotrol.
Similarly in ATAC seq data when I'm comparing two stages how do i define the stage specific peaks? Such as set of peaks are only for A and other for B?
So far what i have done is I get a consensus peaks then from those peaks I get read-count and from those read count I do run deseq2 on those accessible data.
Can I do the same as I do for RNA seq data ,and compare peaks of A vs B using deseq2 and say if a peak is differentially accessible then it can be considered as that condition specific peaks.
Or
There is other way of finding out condition or stage specific peaks ?
"Identification of tissue-specific chromatin accessible regions

We used a strategy described previously based on the Shannon entropy
to compute a tissue specificity index for each peak4,36,37.
Specifically, for each peak, we defined its relative accessibility in
a tissue type i as Ri = Ei/ΣE, where Ei is the RPM value for the peak
in the tissue i, ΣE is the sum of RPM values in all tissues, and N is
the total number of tissues. The entropy score for each peak across
tissues can be defined as H = −1 * sum(Ri * log2Ri) (1 < i < N), where
the value of H ranges between 0 to log2(N). An entropy score close to
zero indicates the accessibility of this peak is highly
tissue-specific, while an entropy score close to log2(N) indicates
that this peak is ubiquitously accessible38. Based on the distribution
of entropy scores, peaks with score less than 3.5 were selected as
tissue-restricted peaks."

The above is one of the way they saying tissue specific peaks from this paper
Any suggestion would be really appreciated

Differential analysis like DESeq2 and limma have been successfully applied to ATAC-seq data. The method you quoted above seem to be too fancy.

i would prefer simple one .  But my question regarding finding peak specific condition or loss or gain of peaks between two condition, can it be addressed through the differential analysis ? I'm just not clear how a loss or gain can be shown in accessibility data through differential analysis

But you mentioned consensus peak set in your question above, so in this way loss or gain can also be defined by differential. condition-specific peaks do not have to be binary loss or gain. Loss or gain may be intuitive for ATAC-seq peaks but the concept itself is confounded by sequencing depth and peak-calling methods, for instance. May be you can do differential analysis first, and choose the peaks that show large fold change, and go back to check if the peak is called in both or one condition.

okay "so in this way loss or gain can also be defined by differential." so is similar what is done for RNA seq lets say A vs B , i did differential analysis so how do i say that some region are specific to A or specific to B  ?


text:
How to tell which channel is the reference channel in two channel (red green) array data?
I'm sorry if this is a stupid question - I have not analysed two channel array data before.
I have a GenePix data file, which does not seem to indicate which channel (red or green) is used as reference channel.
If it possible to tell this from the data itself?
Where raw is my data object, and G is a dataframe of the green channel data and R is the data frame of the red channel data, the ranges are as follows:
#range(raw$G)
[1]    45 21939
>
> range(raw$R)
[1]  39 332

No, you can't get this from the data itself. If there is a control/reference sample on one channel, you need someone to annotate which channel it's on.


text:
freebayes segmentation fault
I have a draft, Illumina-based, genome and bunch of reseq samples. I generated bam files using bowtie2 and marked duplicates using picard.
Now, when I tried to call SNPs using freebayes, I got cryptic Segmentation fault (core dumped) even when I tried to run it in verbose debugging mode (-dd). Do you have an idea where might be a problem?
freebayes -f data/Afus1/genome.fa -b reseq1.sorted.rmdup.bam reseq2.sorted.rmdup.bam -dd
  loading fasta reference data/Afus1/genome.fa
                                                                         Opening 2 BAM format alignment input files

     done
Number of ref seqs: 821466
Number of target regions: 0
no sample list file given, reading sample names from bam file
found 2 samples in BAM file
Segmentation fault (core dumped)
```

Oh you silly sausage!
This obviously means that you are running out of stack memory. Apparently...

freebayes is creating a single string containing the entire vcf header, which gets quite big when the reference has a lot of contigs with long names, causing stack space depletion.

Try increasing the maximum stack size (in this example to 128M) and running your command again.
ulimit -S -s 131072

+1 silly sausage.


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

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

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

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

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

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

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

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

      dev.off();

     cap_ord

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

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

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

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

Many thanks

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

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

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

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

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

Thanks, I have updated the result

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

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

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

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

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

Thanks, I will try to understand these all.


text:
How to programmatically classify a protein according to its genbank feature
Say I found an interesting protein in a genbank file, e.g.:
 gene            273792..275552
                 /locus_tag="ECs_0242"
 CDS             273792..275552
                 /locus_tag="ECs_0242"
                 /codon_start=1
                 /transl_table=11
                 /product="type IV secretion protein Rhs"
                 /protein_id="NP_308269.1"
                 /translation="[...]"

(from https://www.ncbi.nlm.nih.gov/nuccore/NC_002695.2)
I would like to programmatically get a classification of this protein.
The (not new) motivation is that my pipeline finds many proteins, and I want to identify enriched or depleted protein classes (e.g., enrichment for DNA-binding proteins).
By googling for "type IV secretion protein Rhs", I found a UniProt entry which has a link to a QuickGO page, which includes three categories: "transposase activity", "transposition, DNA-mediated" and "DNA binding". From this page I can also get an ancestor chart for each of these categories, e.g.:

Optimally, the code would work quickly (preferably locally), and for a description from a CDS feature (e.g.: "type IV secretion protein Rhs") return the GO categories it belongs to. It would be even better to also have the ancestry of each category in a computer-readable format.
What is the best way to do that?
Is there a specific UniProt subset that fits my need (finding GO categories given a protein description)?
Is there a specific GO subset that fits my need (extracting ancestor GO categories as explained here)?
(I figured that GO categories fit my need of classifying proteins, but if a better option exists, please let me know)

Using the UniProt website API (https://www.uniprot.org/help/api) this could be done with a query like
https://www.uniprot.org/uniprot/?query=name%3A%22type%20IV%20secretion%20protein%20Rhs%22&columns=id%2Centry%20name%2Creviewed%2Cprotein%20names%2Cgenes%2Cgo(biological%20process)%2Cgo(molecular%20function)%2Cgo(cellular%20component)&sort=score
or, in tab-separated format:
https://www.uniprot.org/uniprot/?query=name:%22type%20IV%20secretion%20protein%20Rhs%22&format=tab&limit=10&columns=id,entry%20name,reviewed,protein%20names,genes,go(biological%20process),go(molecular%20function),go(cellular%20component)&sort=score
(Remove the "&limit=10" part from the URL to obtain the full set instead of the first 10 entries.)
You can add colums as described here:
https://www.uniprot.org/help/customize,
https://www.youtube.com/watch?v=ado1r8IDm3U (video, at about 1min25)
or
https://insideuniprot.blogspot.com/2015/03
Note that while UniProt tries to improve protein name annotation, the unreviewed section of UniProtKB (i.e. UniProtKB/TrEMBL) is often far from perfect and its accuracy relies on the information provided by the submitter of the nucleotide entry, unless an automatic annotation rule is available to assign a more suitable name (see https://www.uniprot.org/help/protein_names).
Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.


text:
Why could I be having only NAs for p_val and p_val_adj in Seurat DGE analysis?
I am doing a differential gene expression analysis to find out genes that are differentially expressed in one cell type from the lung in two conditions, control and idipathic pulmonary fibrosis (IPF). In total I have 152538 cells. Can anyone help me provide a hint to understand why only NAs get returned for the pvalues and adjusted pvalues?
> table(seurat$group)

Control     IPF
  64652   87886

> degs <- FindMarkers(seurat,
                            ident.1 = "IPF",
                            ident.2 = "Control")
> head(degs)
         p_val  avg_logFC pct.1 pct.2 p_val_adj
SPP1        NA  1.2721607 0.390 0.125        NA
FN1         NA  1.1768771 0.741 0.409        NA
CRIP1       NA  1.0484934 0.938 0.832        NA
CCL18       NA  0.9830185 0.734 0.469        NA
S100A10     NA  0.5799861 0.965 0.902        NA
VIM         NA  0.5666287 0.992 0.964        NA

are you showing the full function called? it should FindMarkers(seurat,group.by="group",ident.1="IPF",ident.2="Control") unless you set the identity to be group

@StupidWolf Yes, I did set the identity to group. The function works for other cell types in the same data set.


text:
Retrieve overlap widths with Bioconductor
I want to retrieve all overlaps between two sets of genomic intervals (GRanges) using Bioconductor. I am interested in the widths of the resulting overlapping ranges.
In cases with discontinuous ranges, intersect() works, but I am interested in the case where both sets contain contiguous ranges (no gaps).
I am using the findOverlaps function as in the example below:
gr1 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 11, 21, 31), c(10, 20, 30, 40)))
gr2 <- GRanges(seqnames="chr2", ranges=IRanges(c(1, 5, 15), c(4, 14, 30)))
hits <- findOverlaps(gr1, gr2)

This code yields the following intersections in a "Hits" object:
Hits object with 5 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         1           2
  [3]         2           2
  [4]         2           3
  [5]         3           3
  -------
  queryLength: 4 / subjectLength: 3

But I could not find a way to extract the intersections widths from the resulting object.
It this the wrong approach ?

Note: Using intersect(gr1, gr2) in this case yields a single contiguous overlap range spanning the whole chromosome, but I want each overlap to be returned independently, just like findOverlaps does.

Please add a representative output example.

Thanks, I just added it

Whilst intersect is sort of doing what you want, it merges adjacent overlaps together which is not what you want I suppose. There doesn't seem to be dedicated Grange function to do this out of the box, but you could go through all overlapping pairs resulting from findOverlap and calculate it's widths using pintersect:

punion, pintersect, psetdiff and pgap are generic functions that compute the element-wise (aka "parallel") union, intersection, (asymmetric!) difference and gap between each element in x and its corresponding element in y.

"pintersect"(x, y, resolve.empty=c("none", "max.start", "start.x"), ...) # x, y   IRanges objects.

hits <- findOverlaps(gr1, gr2)

width(pintersect(gr1[queryHits(hits)],
                 gr2[subjectHits(hits)]))
[1]  2  6  4  6 10

If I understand correctly, pintersect only computes the intersection between pairs of ranges at the same position. (e.g. 2nd vs 2nd). This implies you have the same number of ranges in gr1 and gr2. I am interested in the all-vs-all intersection. I will update my example to make that clearer, sorry.

Thanks that answers the question :) The only issue is that the R loop is super slow. In my case I have 468,000 intersections and it's still running after 10mins. I have found an alternative which runs in 0.1 sec, but it is imperfect as I have to truncate intervals by 1bp (slightly wrong results):
mywidths <- width(intersect(narrow(gr1, start=2), narrow(gr2, start=2))); comp_df <- data.frame(feat1 = ol@from, feat2 = ol@to, size = mywidths)

That a good point. I've added another possibly quicker solution, hope that works!

That last solution is a great idea, thanks !

If it works for, I'd consolidate the answer since the first iterations are superflous then :)


text:
How to import a large amount of .bed, .gff, .vcf, .paf, .sam files into an SQL database?
Are there best practices to load different bioinformatics file formats such as VCF, BED, GFF, and SAM to SQL databases? I am wondering how people out there do that efficiently.
All of these three formats are tab-separated files, so basically the following should work. I feel weird about it since most people I know don't use MySQL to work with these files.
LOAD DATA LOCAL INFILE 'bed.bed'
INTO TABLE bed-file
FIELDS TERMINATED BY '\t'
LINES TERMINATED BY '\n'
IGNORE 1 ROWS (list of the columns)
SET creation_date  = STR_TO_DATE(@creation_date, '%m/%d/%y');

Did you ever find a satisfactory answer to this question? I'd love to hear more about what sort of solution, if any, you ended up with.

@JessimeKirk not really. I've heard recently about Glow. Maybe it will be helpful to you.

The clodius tool puts tab-delimited data from a BED file into an rtree in a sqlite3 database:
https://github.com/higlass/clodius/blob/c98bb16ade93402fcea3b749d705c52ea165b609/clodius/cli/aggregate.py#L457-L777
Data are stored in this database at multiple "zoom levels", which allows for fast queries at wide genomic extents:
https://github.com/higlass/clodius/blob/develop/clodius/tiles/beddb.py
Aggregation (input) and querying (output) seem efficient in practice, but efficiency may depend on aspects of your specific dataset, i.e., whether you mean efficiency in execution time of INSERT or SELECT statements, or in space efficiency in the overall database size.
Other data can be converted to BED via vcf2bed, sam2bed, gff2bed, etc. and then aggregated. This is not the only way to do things, of course, but hopefully this can suggest directions to explore.

You might also take a look at this genomics-focused sqlite extension: https://github.com/mlin/GenomicSQLite
As described in the repository README:

This SQLite3 loadable extension adds features to the ubiquitous
embedded RDBMS supporting applications in genome bioinformatics:

genomic range indexing for overlap queries & joins
in-SQL utility functions, e.g. reverse-complement DNA, parse "chr1:2,345-6,789"
automatic streaming storage compression (also available standalone)
reading directly from HTTP(S) URLs (also available standalone)
pre-tuned settings for "big data"


text:
How do I infer someone's haplotype of a gene from data on several of the gene's SNPs?
For a a few SNPs from given gene, I know there happen to be several haplotypes that most of the population has. But from a person's data on those SNPs, we have two alleles at each SNP, so it's not trivial to match the SNP data to infer the person's haplotypes. When I search this I find tools that seem to be for haplotype reconstruction and SNP imputation, such as SHAPEIT, but this is not what I'm looking for.
Might someone know of a tool for this inference of known haplotypes from SNP data? It would be helpful even if someone could give me a name for this problem, as I'm not sure of the right term to search.

It depends on the scale of the data (i.e. number of individuals, number of SNPs). Generally this can be done by matching up individual SNP genotypes to haplotype tag SNPs. The specific combination of genotyped SNPs is used to determine the underlying haplotype, based on tag SNP frequencies for each haplotype.
More details of the process of identifying tag SNPs can be found here:
https://www.broadinstitute.org/haploview/tagger

The scale is quite small: only 4 SNPs and ~1k people. Tagger seems to want a local copy of HapMap data to identify tag SNPs, so it's probably overkill. I think it'll be best for me to look at how tagger works and implement my own function.


text:
How to assess how prone genes may be to acquire structural polymorphisms?
I have 5 strains of P.falciparum. Each FASTA file has all its annotated CDSs. After a first pre-processing phase, where I eliminated the strangest sequences (perhaps the longest or shorter ones, which did not start with methionine or which did not have a length multiple of three) I wanted to align the CDS of these strains to evaluate how much every single gene was prone to acquire structural polymorphisms (insertions, SNPs deletions, etc.). A first step I thought could be BLAST to find the orthologs and associate the identifiers of the genes. The problem is that BLASTn has a limit on the number of sequences it can use. How can I go about getting multiple alignments? Also, is there any useful statistic that you can apply that you know of in the literature that can provide a mathematical criterion for obtaining a statistic about it?
thanks a lot

why are you using blastn and not tblastx, for CDS? 2) are you using the blast web portal or a binary on HPC or laptop? I would be surprised if a binary had this limit. 3) have you considered other methods, such as using HMMER to annotate ORFs with Pfam and then matching Pfams to find orthologs? or translating to protein sequences, which would be easier? 4) Have you considered googling for multiple sequence alignment packages?

I Split the comment in two part :I was using the portal. It is the first time that I have to use BLAST massively so I am quite bewildered. Since the web portal has limits in the bases it can use I downloaded the standalone BLAST and was studying how to use it. As for the packages you mention: are you referring to R? Since I have to align 5 proteomes with more than 5000 sequences to find the orthologues among these isolates, if I used R would the computation times be acceptable?

Otherwise what are you referring to? No Pfam I had not considered for now. Aside from the things mentioned above I am adding it to the list of tools to use to get around the problem. Can you link me some reports from which to draw information in such a way as to make it faster to use? thanks a lot

Standalone is likely to work. You will have to use makeblastdb to make a database per proteome. 2) No, these are not in R. tblastx (or blastp if you translate) is a part of the blast suite (http://nebc.nox.ac.uk/bioinformatics/docs/tblastx.html). I suppose that directly searching with BLASTP or whatever will yield best hits, but annotating with Pfam is still something that you will likely have to do anyways. Using hmmsearch (in HMMER) to search against Pfam is probably simplest and fastest, if you choose that route (http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf).

STALN is what I'm currently trying to use now. Indeed, I am creating a DB per proteome. Once I will have studied the material I will post my doubts. Thank you in advance!

I managed to get a file with an output format like in Trace Back Operation (BTOP) like in https://www.ncbi.nlm.nih.gov/books/NBK279682/ . Is there any way to parse a file like that? This is first time I have to deal with such data so I'm sorry for these questions that could result basics.


text:
Why is the UniProt REST API returning multiple results, when I am only providing one ID to be checked?
Can some explain why when I run a REST API query on UniProt for one ID ('Q9UJL9'), the results for three IDs are returned?
The query I'm using is:
https://www.uniprot.org/uniprot/?query=Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab
(I didn't hyperlink it since you need to see the query to see where I'm going wrong).
When I run the query above, the output is:
Entry   Gene ontology (molecular function)
Q9UJL9  DNA-binding transcription factor activity, RNA polymerase II-specific [GO:0000981]; DNA-binding transcription repressor activity, RNA polymerase II-specific [GO:0001227]; metal ion binding [GO:0046872]; RNA polymerase II transcription regulatory region sequence-specific DNA binding [GO:0000977]
Q8WWY3  identical protein binding [GO:0042802]; ribonucleoprotein complex binding [GO:0043021]; RNA binding [GO:0003723]; snRNP binding [GO:0070990]; U4atac snRNA binding [GO:0030622]; U4 snRNA binding [GO:0030621]
Q6A162  structural molecule activity [GO:0005198]

I checked the uniprot page of the query I want (Q9UJL9), the other items are interactors (along with many others) but not synonyms/linked (it seems, e.g my query is a zinc finger, one of the linked records, Q6A162, is a keratin protein).

It makes sense as the search searches for all instances: if you google yourself you are likely to find photos of your colleagues, because their appear in the same page as your name, same thing here.
If you already know your protein an alternative is to get https://www.uniprot.org/uniprot/Q9UJL9.txt and grep out the DR   GO;  entries, but that is as convoluted.

it is because you do a text search and the "word" Q9UJL9 appears in multiple entries. If you want an id search, specifically say so and for that the query syntax is id:Q9UJL9. You will want to follow a redirect, and there are some rare cases where that might still lead to more than one entry if the id was made secondary and attached to multiple entries. e.g. P00001.

This is perfect thanks a mil I totally get it now.


text:
Quast duplication ratio and mismatches percent
When analyzing Quast results it seems that it doesn't calculate mismatches and indels in a useful way if the "Duplication ratio" is over 1.
For example, that's what I get for an assembly with two contigs covering almost the full genome

And that's what I get if I remove one of the contigs:

See how the mismatches and indels per 100kbp are almost double in the first case but that does not correspond to the reality.
Any better way to get these metrics? Would it be safe to divide error rates by the duplication ratio?

I took a quick look at the quast code that I believe is responsible here:
def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
    indels_info = IndelsInfo()
    genome_mapping = {}
    for chr_name, chr_len in reference_chromosomes.items():
        genome_mapping[chr_name] = [0] * (chr_len + 1)
    with open(used_snps_fpath, 'w') as used_snps_f:
        for chr_name, aligns in ref_aligns.items():
            for align in aligns:
                ref_pos, ctg_pos = align.s1, align.s2
                strand_direction = 1 if align.s2 < align.e2 else -1
                for op in parse_cs_tag(align.cigar):
                    if op.startswith(':'):
                        n_bases = int(op[1:])
                    else:
                        n_bases = len(op) - 1
                    if op.startswith('*'):
                        ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                        if ctg_nucl != 'N' and ref_nucl != 'N':
                            # MISMATCHES INCREMENTED HERE
                            indels_info.mismatches += 1
                            if qconfig.show_snps:
                                used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += 1
                        ctg_pos += 1 * strand_direction
                    elif op.startswith('+'):
                        # INDELS INCREMENTED HERE
                        indels_info.indels_list.append(n_bases)
                        indels_info.insertions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = '.', op[1:].upper()
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ctg_pos += n_bases * strand_direction
                    elif op.startswith('-'):
                        indels_info.indels_list.append(n_bases)
                        indels_info.deletions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = op[1:].upper(), '.'
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += n_bases
                    else:
                        ref_pos += n_bases
                        ctg_pos += n_bases * strand_direction
                if align.s1 < align.e1:
                    for pos in range(align.s1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
                else:
                    for pos in range(align.s1, len(genome_mapping[align.ref])):
                        genome_mapping[align.ref][pos] = 1
                    for pos in range(1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
            for i in ns_by_chromosomes[align.ref]:
                genome_mapping[align.ref][i] = 0

    covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
    return covered_ref_bases, indels_info

This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.
This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).
From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast used_snps file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.


text:
How does one convert a log-likelihood substitution matrix into a probability matrix?
I am trying to calculate the substitution probabilities for amino acids (i.e. the probability of one amino acid given another). In theory, you should be able to derive them from a log-likelihood matrix (also called log-odds), such as Blosum62.
Does anyone know how to calculate probabilities from log-likelihoods or have a good modern source for a substitution probability matrix?

If I understand what you want I think you just need to align the sequences and then make some probability calculation (which will likely make no sense biologically but will work for you). If the sequences have rearrangements as I understand you said somewhere then it's more complicated than just a multiple-alignment and you need a tool that finds local alignments. Am I following you?

I mean, if I understand well that's not going to be possible with any probabilistic method. You need to align first somehow. I guess something similar to what QUAST does. No idea if there is a tool for that. I said...

Yes, the sequences are aligned first before comparing the amino acid probability of each position between the two groups. I was just wondering how to convert a Blosum62 log-likelihood score of "2" (for example) to a probability, because I can then use the set of probabilities to calculate a combined probability at every position for a set of sequences.

Note that you are using imprecise language. You don't specify probability of what.  Anyway, I don't think I can help with that part.

UPDATE2
@chubs wants to compare different clades. This is possible within recent releases of ConSurf:

It is additionally quite possible to use ConSurf's default outputs to get per-position conservation scores that can be used for post-processing in any fashion, as well as prefab scripts that can be used with PyMol to make custom visualization:

For all cases, ConSurf creates the following outputs:

The sequence and MSA colored by ConSurf conservation scores.

A text file that summarizes for each position the normalized score calculated, the assigned grade (1-through-9), the reliability estimation (for the Bayesian method) and the amino acids/nucleotides observed in the respective MSA column.

The sequences selected for the MSA and the MSA constructed (unless those files were uploaded by the user).

A file with the frequency of each amino acid/nucleotide observed in each column of the MSA.

The evolutionary tree, which was calculated by the server or uploaded by the user, together with the MSA are shown using the WASABI platform (Veidenberg et al., 2015).  Moreover, using WASABI, a user can select a subtree containing a fraction of the homologous sequences and conduct a follow-up ConSurf analysis with these selected sequences. To refine a ConSurf analysis to a selected subtree, a user can choose any internal node on the phylogenetic tree and open a WASABI menu by using a right mouse click. Selecting the option 'run ConSurf on subtree' will issue a new window with the new ConSurf run for the selected subtree sequences (see example in Figure 3).

UPDATE
It seems that what @chubs really wants to do is map protein conservation onto a structure. There are various tools for doing this (here, here). ConSurf is the older and more established tool, using evolutionary conservation to rule out trivial apparent associations due to phylogenetic non-independence.
here is the  somewhat related 3DPatch workflow, which uses profile HMM information as a measure of conservation:

Here is the output that they get:

Likewise ConSurf gets similar looking outputs:

ORIGINAL ANSWER
I think that more information is needed here to answer this question properly (what is your data? what is your ultimate goal? what are the things you've already tried?), but we can say a few things right away.
I would strongly suggest looking into substitution models, which are literally substitution probability matrices. Note that these are probabilities per unit time, which is to say that they are calibrated to e.g. evolutionary rates along a specific phylogenetic tree. I would therefore strongly suggest incorporating a phylogenetic model of some kind.
If you have e.g. a multiple sequence alignment, you can estimate a phylogenetic tree by one of the many tree estimation methods, which will give you a substitution probability matrix "for free" as this estimation is taken care of in the process of tree estimation. These will again be substitution probabilities per unit time; you will honestly not be able to get rid of the time-dependence here, as that's how evolution works.
I am a little behind the times, but some tools you might look into for tree estimation are PhyML and FastTree. If you only have sequences and haven't aligned them yet, you can use a tool like MUSCLE, ClustalO, or MAFFT.
BLOSUM and related models are derived from specific slices of protein sequence identity (evolutionary history) calibrated to certain scopes, it's unlikely that they are appropriate to your case, unless your data fits inside those rather narrow slices.

I have a set of aligned sequences for which I want to generate a probability distribution of amino acids at each position. For example, at position 1, if all amino acids in the alignment are alanines, that position would hold a probability distribution for just alanine, which would be high for itself (alanine) and lower for other amino acids (just like the blosum62 row for alanine). If a position has different kinds of amino acids, its probability distribution would a mix, accordingly.

My ultimate goal is to have two seperate sets of sequences that can each be distilled to a probability distribution at every sequence position. I can then compare the two distributions (one from each set) to determine how similar that position is between the two groups.
Sorry if this is confusing, but it's the best way I can think of for comparing two sets of alignments (e.g. mammalian sequences versus their arthropod counterparts).

@chubs is there a reason that you do not just use the sample distribution for each column of the alignment (for each taxonomic slice of alignment rows, if needed)? This is the normal methodology for e.g. generating logo plots, which I would recommend looking into: https://en.wikipedia.org/wiki/Sequence_logo If you are feeling very Bayesian you could pseudocount or try to incorporate a prior distribution based on the background frequency of amino acids in all positions of all the sequences, though that is getting a little over-elaborate IMO.

@chubs I'll note that methods such as PAML have specifically implemented phylogenetic modeling approaches for detecting amino acid patterns of the type that you seem to be interested in between clades, with the branch model allowing $\omega$ variation among branches: http://abacus.gene.ucl.ac.uk/software/paml.html

The reason I'm going about it this way is because I want to map the similarity of a given amino acid between two groups onto a three-dimensional atomic model of the protein. I've tried using sequence identity, but the results aren't quite what I want. In other words, one group with all glutamates at position 1 versus another group with all aspartates at the same position will give a 0% similarity when in fact this conservative substitution should give a higher score. I thought that by using probability distributions, I could calculate a more accurate similarity score to map onto the structure.

see edit to my answer.

Consurf is great but doesn't compare two groups to each other, right? It just calculates the conservation of a set of sequences. I actually want to compare two groups to each other (e.g. multiple mammalian sequences in group 1 to multiple arthropod sequences in group 2). That's why I'm trying to calculate the probability distribution at each position: to then calculate a score based on how similar they are between the two groups.

@chubs then why not use consurf on each clade independently, and compare? Alternately, it seems like it should be possible to use any of various tools to export values for each position in each clade of interest, assuming that they are roughly structurally similar. Then you can do whatever you want with the conservation estimates, including mapping them onto a structural model. PAML would also work for this of course, with the advantage that it is explicitly modeling divergence between your two taxa. But you could also fit profile HMMs for each clade and compare, for the conservation task.

@chubs see incoming update- you can run consurf on subtrees.

If I use consurf on each clade (or on a subtree) what I will get is how conserved an amino acid is within each clade (e.g. position 1 of clade 1 is 10% conserved and position 1 of clade 2 is 30% conserved). This isn't what I'm looking for. I want to know how similar the amino acids are between the two clades (e.g. position 1 is 60% similar between clades 1 and 2). The two clades don't always have the same number of sequences, so I thought I would calculate the probability distribution at each position to compare instead of running Consurf on the combined sequences. Sorry if it's unclear.

@chubs I see. I continue to think that this is a weird way to do it, and you should use a phylogenetically informed approach such as PAML, which will give you a more exact estimate of what is going on at each position in each clade. It seems like that isn't what you want, though. It's quite heuristic, but if you like, you should be able to do a  summation across the amino acid frequency distribution at each site with BLOSUM (or whatever) as a transformation. Do it as a matrix multiplication of frequency vector 1 * BLOSUM * frequency vector 2. That is a per-site similarity measure, sorta.

@MaximilianPress But I guess he has to align first...If I understand what he says it's not always ordered If that's the case that will be the hardest part.

But I agree that the phylogenetic tree way would be the ideal way (despite that also requires a previous alignment)

@juanjo75es based on what they says ("I have a set of aligned sequences for which I want to generate a probability distribution of amino acids at each position.", first comment on this answer), I think they already have the alignments? Maybe I am misunderstanding something though.

@MaximilianPress "The two clades don't always have the same number of sequences" "one group with all glutamates at position 1 versus another group with all aspartates at the same position will give a 0% similarity" That doesn't seem to mean that it's already aligned.Or at least not well-aligned. No idea. If he says it is then let's accept that.


text:
installing multiple bioconductor packages at once
I was wondering if there is more elegant way of installing and loading multiple packages in Bioconductor similar to pacman with CRAN packages.
I tried:

# install and load  the package  manager
 if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# install and load bioconductor packages
bio_pkgs = c("BiocManager", "ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "cBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

 if (!requireNamespace(bio_pkgs, quietly = TRUE))
    BiocManager::install(bio_pkgs)

library(easypackages)
libraries("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData", "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

You can do something like this:
bio_pkgs <- c("ggbio", "biomaRt", "EnsDb.Hsapiens.v75", "renvcBioPortalData",
          "GenomicRanges", "GenomicFeatures", "TxDb.Hsapiens.UCSC.hg19.knownGene", "Homo.sapiens")

# install:
BiocManager::install(bio_pkgs)

# load all at once
invisible(lapply(pkgs, function(x) library(x, character.only=TRUE)))

Be aware though that if you version-control your packages with something like renv then it is encouraged to explicitely load every package as suggested here so renv can properly recognize the packages to snapshot.

Thank you ATpoint, this is way more elegant :)


text:
Find pattern that is present twice and allow <=2 mismatches on each
I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with <= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-

TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA

Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.

Do you have any additional assumptions you can make? For example, are the pairs of barcodes within a specific distance from each other  or the reads ends? I ask because with an edit distance of two (including indels) the chances of random matches will start to increase

This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.
I just copy pasted your examples seqs into files to read (called them seqs and code).
#need blast installed
import subprocess
import pandas as pd

path = '/Users/ur/dir/'
#make a blast db
with open (path + 'seqs.txt', 'r') as fin:
    with open (path + 'seqs.fa', 'w') as fout:
        for i, line in enumerate(fin):
            fout.write('>'+str(i)+'\n')
            fout.write(line)

subprocess.run([
    'makeblastdb',
    '-in', path + 'seqs.fa',
    '-dbtype', 'nucl'
])

#use blast
subprocess.run([
    'blastn',
    '-task', 'blastn-short',
    '-query', path + 'code.txt',
    '-db', path + 'seqs.fa',
    '-out', path + 'out.txt',
    '-outfmt', '6'
])

#grab seqs with 2 hits (ie ~2x len(query)
df = pd.read_csv(
    path + 'out.txt',
    sep='\t',
    header=None,
    names=['query', 'db_name', 'percent','len'],
    usecols = [0,1,2,3]
)
df = df.groupby('db_name').agg('sum')
print(*list(df.query('len > 27').index))

I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.


text:
how to easily make complex Euler/Venn diagrams in R?
From my microarray data I have selected three subsets of probes. I would like to display these subsets as a Euler diagram to visualize the overlaps, but am struggling with inputting the combinations into Eulerr. All three subsets need to be nested within the full set of probes; one subset is fully contained within another, and the third overlaps with the other two. Is there a package into which I could feed the vectors of probe names and let it determine the overlaps for me? Many thanks in advance!

Since you are in R, consider using the UpsetR library.  It doesn't make Venn diagrams but it helps to visualize overlaps between any number of groups.
http://gehlenborglab.org/research/projects/upsetr/#:~:text=UpSetR%20is%20an%20R%20package,based%20on%20groupings%20and%20queries.

You might consider different ways of using the venneuler package. I'm not super thrilled with its usability (or with the accuracy of the overlaps), but you can get it to make custom Venn diagrams.
Here is a thread with some suggestions along these lines.
Below is some sample R code that I slightly modified from an old project that I used to make a 3-category Venn diagram.
I didn't run this to make sure that it works because of an irritating Java dependency issue on my OSX laptop that I don't feel like fixing right now, but as long as you have a JDK on your machine and can install the packages it should work.
# note irritating rJava dependency of venneuler
library(plotrix)
library(venneuler)

#A-B sharing: 200414
#B-C sharing: 17561
#A-C sharing: 16764
#All share 14201

# get coordinates for circles- replace as needed
# this function is supposed to plot stuff, but doesn't actually
ven = venneuler(c(C=48147-17561-16764+14201,
                  B=347181-17561-200414+14201,
                  A=325971-200414-16764+14201,
                  "C&B"=17561, "B&A"=200414, "C&A"=16764, "A&B&C"=14201))

cats = c('A','B','C')

# note- this following only shows up in the .pdf printer, not in the X11 utility plotter
pdf('venn.pdf')
plot(c(-.5,1.5), c(-.5,1.5), type='n', xaxt='n', yaxt='n', bty='n', xlab='', ylab='')
draw.circle(ven$center[1,1], ven$center[1,2], ven$diameters[1],
        col = hsv(1,1,1,.5))
draw.circle(ven$center[2,1], ven$center[2,2], ven$diameters[2],
            col = hsv(.5,1,1,.5))
draw.circle(ven$center[3,1], ven$center[3,2], ven$diameters[3],
            col = hsv(.1,1,1,.5))

text(ven$center[1,1]+.065, ven$center[1,2], labels=cats[3], cex=1.3)
text(ven$center[2,1], ven$center[2,2]-.2, labels=cats[2], cex=1.3)
text(ven$center[3,1], ven$center[3,2]+.2, labels=cats[1], cex=1.3)

dev.off()


text:
How do you determine the strandedness of a RNAseq dataset?
How do you verify if a RNAseq dataset is unstranded, stranded or revesely stranded from a fastq file?

I was looking for a method that didn't involve quasi-mapping but I think this will solve my problem. Thank you!

Unless the metadata is recorded somewhere you won’t get around mapping: the FASTQs themselves don’t record this information anywhere.


text:
How can I extract information from .sam files?
I have 10 .sam files after my bowtie2 alignment on ten single-pair sequences. I would like to build a graph based on that output data, however I don't see any other solution but just writing the alignment percentage manually after each iteration.
Is there a tool that can help me to build a .csv file so I can manipulate the output data more easily?

You should never have SAM files: they’re not useful, not even as intermediate files! You should always directly generate BAM files. They’re smaller, more efficient, and allow more features (e.g. indexing for even faster access and statistics). For example, you would usually call bowtie2 as follows: bowtie2 ‹args› | samtools sort --output-fmt-option level=0 | samtools view -b -o sorted.bam (this immediately sorts the BAM file, which requires a lot of RAM; if this isn’t suitable for you, omit the intermediate step; but you normally want to have sorted BAMs).

Use samtools flagstat with option -O tsv:

Using -O tsv selects a tab-separated values format that can easily be imported into spreadsheet software. In this format the first column contains the values for QC-passed reads, the second column has the values for QC-failed reads and the third contains the category names.

SEE ALSO:
What Does Samtools Flagstat Results Mean?
INSTALLATION:
Some old versions of samtools do not have option -O tsv. To install a recent versions of samtools, use either:
conda (macOS, Linux) (but see also: Tutorial: installing samtools with conda - overcoming two common errors):
conda create --name samtools samtools

or brew (macOS):
brew install samtools

Thank you, it worked however just in a "default" combination, e.g. samtools flagstat input.sam           If I try  samtools flagstat -O tsv input.sam, it throws me an error  invalid option -- 'O'

@DmitriiTrubetskoy See updated answer.

thanks, however my  samtools is already the newest version (1.7-1). (info from terminal). Just to be sure, the right syntax will be: samtools flagstat -O tsv input.sam ?

@DmitriiTrubetskoy Exactly, and note that 1.7* is older than 1.11, which is the most recent version available through both conda and brew (see links in my answer). I also confirmed that this option is present in samtools 1.10 that I installed via brew. So try updating/installing the new samtools version if you really need this option.

Thanks a lot. I just wonder how can I roll back to 1.10? Is there a way to do it?

@DmitriiTrubetskoy Don't roll back, because 1.7 is fairly old. I am suggesting updating from 1.7 to the most recent version, currently 1.11. Use conda or brew. To install conda (for mac or linux), follow directions here for miniconda: https://docs.conda.io/projects/conda/en/latest/user-guide/install/ . To install brew (for mac), follow the directions here: https://brew.sh/ . Then use the commands in my answer to install samtools (via conda or brew). This will install by default the most recent version. You are usually best off with the most recent version. Good luck!

Thank you. Now it worked, something happened...however I suppose I must have created some new file that would contain some output information - what kind of file that should be? I cannot see find anything in the same folder where input.sam file was.

@DmitriiTrubetskoy It writes into STDOUT a tsv file, 16 rows x 3 cols.

Great, now everything's ok. Much appreciated the tips and for being patient.


text:
What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)?
I bought a "Whole Genome Sequence" kit, which gave me the following VCF files (this service did not include any analysis, just the raw data):
cnv.vcf
indel.vcf
snp.vcf
sv.vcf
My understanding is that these files show the differences of my DNA relative to a reference human genome (hg19 in my case)
What tools and resources can I use to find out which known alleles/variants of genes I have?
An example: Blood type
According to Wikipedia, there are 6 common alleles of the ABO gene, and which alleles I have determine my blood type.
Looking at my VCF results, I have these SNPs and INDELs within that gene:
<removed>

Given this, how would I figure out my blood type?
I understand that my data may not be 100% accurate, and that many phenotypes are the result of more than one gene, environmental factors, or are not well studied at all, but this should work for some simple things ie blood type, COMT, BHLHE41, CYP1A2, etc

For the blood locus. Try and find online the genomic position (i.e. the number like 133273540 in the POS field of the .vcf) of the ABO blood locus and look it up on your dataset. Remember to make sure you find it on the gr37/hg19 build.

Wait, did you got yourself yourself sequenced?? I thought "my DNA" means "my DNA sample" not "my own DNA"!

And do I get right, you know your variants and you would like to know something about them, right?

Please don't just share your own genetic information online, especially if you are not fully understanding what you are looking at. It is not something to be posted everywhere, although some loci might be relatively harmless you shouldn't take this lightly.

Please ping me if you decide you want this censored. You should be very aware of the long-term consequences for you and your current and future relatives of sharing your genetic information on the internet in a post.

@DevonRyan thanks, I've removed it. I didn't think a small snippet would be a problem

@WouterDeCoster also thanks for the warning

If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.
What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.
There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.


text:
Error in if (is.na(sra_acc$run[i])) { : argument is of length zero
I am trying to write a code for importing sra files as fastq files in Rstudio (for obtaining and aligning RNA-seq reads from public repositories). I keep getting "Error in if (is.na(sra_acc$run[i])) { : argument is of length zero" right after I run this code. Any suggestions on what might be going wrong? My code is as follows:
library(SRAdb)
sqlfile <-'SRAmetadb.sqlite'
if(!file.exists('SRAmetadb.sqlite')) sqlfile <<- getSRAdbFile()
sra_con <- dbConnect(SQLite(),sqlfile)
sraInf <- getSRAinfo("SRP045534",sra_con, sraType="sra")

Error in if (is.na(sra_acc$run[i])) { : argument is of length zero

Have you confirmed whether file.exists('SRAmetadb.sqlite') is TRUE or FALSE? Or generally, what the value of any of your variables is? I am not familiar with this library but having sqlfile be a string that is overwritten by the return value of getSRAdbFile() (don't know what that is) seems like it might cause an issue. Also, are you intending to use the superassignment operator (<<-)? I believe that you are using it in a non-standard way: https://stat.ethz.ch/pipermail/r-help/2011-April/275905.html

After checking with file.exists(), I would suggest you try the same command with another SRA run. It may happen that some runs are not accessible for reasons I do not know. I always have problems in accessing sra formatted files and use a workaround to download fastq.
Here is my not official workaround:
runlist<-sraConvert("SRP045534",sra_con=sra_con)
runlist[1:2,]

study submission    sample experiment        run
1 SRP045534  SRA179382 SRS685546  SRX681992 SRR1552451
2 SRP045534  SRA179382 SRS685544  SRX681990 SRR1552449

for (dloadrun in runlist$run)
{
    getSRAfile(dloadrun,sra_con=sra_con,fileType="fastq")
}

This should iteratively download all the runs associated to the project you are studying.
It is possible that you can do the same with sra type file, but I have never been able to use them with sraDB, and never investigated the issue.

I am also having this issue- my SQL file comes back as TRUE when I use file.exists but then in Fabio's workaround, the table comes out blank. Same error as you if I use getSRAfile.


text:
Building small molecules de novo
Can anyone recommend a protocol for building the 3D structure of a small molecule de novo, given only a 2D rendering? So far my idea has been to build the molecule using the fragment builder in PyMol, then make conformers using Frog2. But I have no idea how realistic the starting structure is from PyMol and how exhaustively the conformers are sampled in Frog2.
The downstream application of this would be to dock the ligand into a protein binding pocket. If anyone has a more rigorous protocol for making sure all the conformational space is fully sampled, it would be greatly appreciated!

Most universities have a ChemDraw site lincence, so you may want to check that for drawing a non–3D-embedded molecule. Alternatively there are many other programs. PyMOL is not really a good choice. Also there are online drawing tools everywhere, where you can download a SMILES (see below).
Chemical formula is the wrong way to go, unless you are "enumerating".
The chemical formula of glucose is C6H12O6. This is the same for all hexoses. What you want is a SMILES. This string gives the connectivity of the atoms. For D-glucose it's OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O, the @ symbols are chiral centres.
In Cheminformatics making a list of all possible compound with n-carbons is a big deal, this is called enumeration and you likely don't have to reinvent wheel. Furthermore, there are datasets that are filtered by the Lipinski rule.

GDB datasets (if you don't care if you cannot buy it)
Zinc-DB  (if you care if you can buy it)
Enamine Real (if you care if you can buy it reasonably)

Thank you for the advice. I should clarify, the chemical formula is not what I am starting with, I do have the 2D structures rendered in ChemDraw with the absolute stereochemistry. My problem is really, how to best sample torsional angles at all rotatable bonds to get a reasonable 3D structure?

Rendered as in exported as a raster image (JPG) or 2D draw (as opposed to 3D embedded)? If the latter, in ChemDraw, right clicking on the molecules allows one to export as a SMILES string, which is can go into one's notes or whatever. For downstream applications, exporting as a mol file is the best bet.
To generate conformers, open-babel is an excellent option and has some ways to generate good diversity of conformers —including ones with unhappy torsions if one wanted. If the ligand is a natural metabolite (e.g. ATP) the biggest caveat docking is getting the protonation right.

I should say that compchemists generally used the Python module RDKit —but that requires knowledge of Python.

I have them as a JPG image, that's why I used PyMol to build them and get the stereochemistry correct (it's a bit more intuitive for me to see them in 3D rather than to write the SMILES string by hand). I then exported them as a mol file. I have tried both OpenBabel and RDKit, both seem good, I'm just not sure how extensive the conformational search is


text:
Splitting a GenBank file into smaller files
I have a genbank file that I have been trying to run in HMMER using the following command:
hmmscan --domtblout NC_000962.3.final.gbk.out -E 0.000001 prot_profile.hmm NC_000962.3.final.gbk

However, since the file is too big for HMMER, I get the following error:
Fatal exception (source file p7_pipeline.c, line 697):
Target sequence length > 100K, over comparison pipeline limit.
(Did you mean to use nhmmer/nhmmscan?)
Aborted (core dumped)

I need to split this file into smaller files for the HMMER run. The problem is, I would need the ORIGIN field (containing the nucleotide sequence) at the bottom of the file split into relevant chunks into the subfile as well. This is where I am getting stuck.
I have already looked at GenBank Parser and seqretsplit to no avail.
I would greatly appreciate some help with this. Thank you.

Can do it with biopython

Could you please expand on your answer, Liam?

You can use NCBI EntrezDirect to download data in GenBank format for a specific region of the sequence as follows:
efetch -db nuccore -id NC_000962.3 -format gb -seq_start 1 -seq_stop 99000 > seq.gb

The genome fragments from NCBI do not fit my purpose. The input file is obtained from another program and I need to split that file as mentioned in the question.

I don't have time to give you a polished fully functional version but this will get you started. The next step will be a bit harder as you will need to capture the start of the first CDS in the chunk and end of the last CDS instead of just hardcoding as I have here.
fout = open('/path/first3.gbk','w')
for record in SeqIO.parse('/path/NC_000962.3.final.gbk','genbank'):
    for i, feature in enumerate(record.features):
        new_record = record
        new_record.seq = new_record.seq[0:4437]
        new_record.features = new_record.features[:3]
        if i > 2:
            SeqIO.write(new_record, fout, 'genbank')
            fout.close()
            break


text:
Clustering RNAseq fold-change data
I have a dataset that looks something like this

Treatment1
Treatment 2
Control

Sample1
3.23.
0.87.
1

Sample2
1.71.
1.
1

Sample3
2.88.
5.65
1

Sample4
0.77.
1.34
1

The numbers describe a fold change with respect to the control (where the control is always set to 1).
Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?
By extension, does it also make sense to compute something like Euclidean distances (e.g. fold change vector of Sample1 vs. Sample2) for clustering fold change data? I think clustering on counts/TPM makes more sense, but wasn't sure about doing it on fold changes.
Something tells me that it's not possible because the fold changes are relative values as opposed to measured numbers. I've opted for looking at correlations (e.g. Spearman) for now though; are there alternative metrics/strategies to consider?
If there are some good papers talking about clustering and manipulating fold change data I'd be equally grateful, thanks!

Your control is 1? How did that happen? Would it be more straightforward to cluster the samples/genes directly on the expression data (often visualized using heatmaps). e.g.: https://www.ebi.ac.uk/training-beta/online/courses/functional-genomics-ii-common-technologies-and-data-analysis-methods/biological-interpretation-of-gene-expression-data-2/

Ah, it's been given to me where the control is deliberately set to 1 (as it's differential expression with respect to the control?)

Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?

This is problematic, if you want to compare treatments then compute the actual fold change and its significance. You can't really base anything off of the fold-change (log or otherwise) alone.

does it also make sense to compute something like Euclidean distances for clustering fold change data

It's not unheard of to do this. The problem is you lose all information about how reliable those fold-change estimates are. This is something maintained when you cluster sample-level information (e.g., TPM or robust log transformed data). Not that using variance-stabilized data for clustering will work better, since it tamps down on how much highly-expressed genes influence your clustering (assuming you don't want that to occur).
Have a look at the DESeq2 vignette for some nice examples of what good practices are.

Great, thanks Devon! Having come from a different branch of bioinformatics, it's not exactly clear what to do next once you have those DEGs. I agree that clustering on DEGs is, well, doable, but I can't tell if it's the "right" thing to do as I feel you're doing this on relative data (i.e. ratios) as opposed to an actual number.
e.g. if we applied clustering by Euclidean distance, what does it mean when you have to subtract fold changes?
That was my big debate, really... I'll have a look at that vignette!!