text:
Finding smaller sequences from within larger sequences
I am currently working with fastq files which have hundreds of thousands of lines of text. However, not all of them are sequences I am interested in. My sequences are in one line and have a fixed length of 272 characters.
I possess a list of smaller sequences, around 2 million text strings in the following formats (which are fragments from the various fastq files I have):
Format 1
ATCGATCG
ATCGATCGATCG
ACGATCGTGTTACG
ATTTACGTACGTA
AATCG
ACGATACGATACG

Or I also have them in
Format 2
ATCGATCGATACGAT|ATACGTGTTACGAT|ATCGATACG|ATACGATGA

(I just wrote random sequences just to give an idea of how the pattern.txt files look like, I know Format 1 and 2 are not the same in this post, however, the files are in the same order - only the syntax is changed.).
What I've been trying to do is use grep in a Linux environment and use the pattern.txt files to search for the bigger sequences in my fastq files where the pattern is present, and then output it into a .txt file to further work on.
I've used this command so far (with the Format 2 of patterns):
grep -E -f patterns.txt target.fastq > output.txt

Now, this works.
However, I have to limit the number of patterns in my key file and the number of fastq files I can search at once, otherwise I get "grep: memory exhausted" error.
I also want to print out all the sequences up to and including my pattern, but not past it (I don't care what's after my pattern sequence).
And for that I have successfully used the following:
grep -oE "\S+singlepattern" target.fastq > output.txt

However, I can only manage to get this to work if I manually write in the command the pattern sequence I want to look for (thus I can only use a couple). I can't seem to get it to work with the file list.
So I thought of using a while loop.
This is what I've written so far in the terminal (it's one continuous line, but I formatted it here for simplicity on the eyes):
while read line; do
grep $line target.fastq > output.txt;
done < patterns.txt

However, this does absolutely nothing and just spits out a blank output.txt
I have tried changing the syntax, with quotations on arguments, no quotations, I tried the various -E -f -o -F, egrep, etc but it's all the same. No error message, just blank file.
I've also tried implementing the \S to the loop as follows:
while read line; do
grep -oE "\S+$line" target.fastq > output.txt;
done < patterns.txt

Yet, the result is the same.
Any suggestions on how I can fix this? Or even other methods to get the cut-down sequences (up to and including my pattern) from the fastq files into an output.txt file?
Thanks in advance!

Could you please add some more context to explain specifically why you want to use grep? There are lots of sequence mapping programs that are specifically designed to quickly look for DNA sequences within other DNA sequences (both with and without errors), and there's a good chance your problem could be solved by using one of them.

Searching 2 million patters (m) over hundred of thousands line file(s) (n) is costly O(mn) and can be memory-wise too. I would first de-duplicate my list of patterns, if luckily they are highly duplicated, this can improve performance greatly. grep loads one line at a time in memory, if you use the format 2 as input for patterns, you are effectively loading the whole thing into memory at once and thus the memory issue can present. I would therefore use Format 1, or run it in a server with high memory available

Hello! Thank you for your advice. I produced the pattern file from R Studio and initially they were 2 million individual lines, however, I managed to output a new pattern.txt file with just unique patterns and now it's down to 700k. I do get the "grep: memory exhausted" nonetheless, thus why I'm trying to create a loop :/

@gringer I was recommended grep by a colleague, as well as the fact that most of my desired files are on a linux computer cluster. My laptop is not capable of running the required programming and I have zero knowledge of other software/programs as this is my first time doing dry lab work. Any recommendations would be appreciated and I can have a look nonetheless and see if I can use them!

@JeanLuc I was not able to reproduce the "grep: memory exhausted" error message. I am running in a server but I think PC should do fine too.
I put together a file with repetitions of your Format 1 patters so that I have 1M total patterns, one per line (patterns.txt). I added the regex you need as a prefix to each pattern (eg. \S+ATCGATCG).
Then I took a bzipped fastq file and did 'zgrep -o -f patterns.txt file.fastq.gz'. It works. Slow as hell but that is expected given the complexity of the task (each time you search one pattern across all lines in the fastq file, and again ..)

@JeanLuc for faster lookups, you need not grep but some more complex tool that builds a DAWG or an index over the sequence lines of the fastq file first,

Thank you all for your help. I managed to fix it in the end. The error was with how my pattern.txt file was formatted.
I initially exported the .txt from R Studio by specifying the EOL as \r\n. However, I think this was the issue as then from my laptop to the computer cluster the files seemed to have a windows EOL.
Thus on the linux server I did
dos2unix patterns.txt

And since then all of the suggested loops, both while read and the for chunk work beautifully and I get a nice \S+$line exported into the desired output.txt.
Thanks again all for your advice and help!

For future reference, "\n" is a better line ending to use when creating files from within R/RStudio, especially on a Linux system.

Yep, that's why I asked if you could have Windows style line endings. Glad you sorted it out!

The reason your loop is giving a blank file is because you are using > which will overwrite whatever was previously in the file. This means that your output.txt will only have the result of the last grep and if that last grep has no results, you'll get an empty file.
Now, you're on the right track here. You could use a loop, although it will be horribly, horribly slow (by the way, you can use this exact format in the terminal as well, it doesn't need to be on the same line):
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt

while read line; do
    ## Use >> instead of > so you append to the file
    grep -oE "\S+$line" target.fastq >> output.txt;
done < patterns.txt

Your attempt with the -f option is the most reasonable. You might get it to work if you tell grep to treat the input as fixed strings and not regular expressions (-F). Also, you don't need the -E (and that might be causing issues), so try this in case it works:
grep -Ff patterns.txt target.fastq > output.txt

If that doesn't work, just split your patterns.txt file into smaller files and then use grep in a loop with each of them:
split -l300000 patterns.txt pattern_chunk

That will split patterns.txt into multiple files of 3000000 lines each, named pattern_chunkaa, pattern_chunkab etc. For example, this is the result of running that command on a file with 5 million lines:
$ wc -l patterns.txt
5000000 patterns.txt
$ split -l300000 patterns.txt pattern_chunk
$ ls
pattern_chunkaa  pattern_chunkad  pattern_chunkag  pattern_chunkaj  pattern_chunkam  pattern_chunkap
pattern_chunkab  pattern_chunkae  pattern_chunkah  pattern_chunkak  pattern_chunkan  pattern_chunkaq
pattern_chunkac  pattern_chunkaf  pattern_chunkai  pattern_chunkal  pattern_chunkao  patterns.txt

So, you can now use each of those as the patterns:
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt
for chunk in pattern_chunk*; do
    grep -f "$chunk" target.fastq >> output.txt;
done

This doesn't help you get the sequence only up until the matched pattern. For that, you can add the \S+ to the beginning of each pattern in patterns.txt (this assumes you are using the file which has one pattern per line) and then split the file and loop over the chunks again:
## add the \S+
sed -i 's/^/\\S+/' patterns.txt
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt

for chunk in pattern_chunk*; do
    grep -Ef "$chunk" target.fastq >> output.txt;
done

Note how I keep using > output.txt. This is not actually necessary, you can use command > file even if file doesn't exist. However, it is useful here to make sure that the output.txt file doesn't have the result of a previous run since we are now appending (>>) instead of overwriting (>) its contents.

Hello! Thank you for the answer but it's still not working. The grep commands all work, but as soon as I write them within the while or for loops they don't work. To be sure, I also added do echo $line or $chunk (whether it was while or for loop), and it shows that it's running through  all the files. However, the output files are always blank EDIT: I wrote the echo both before the grep and after, and it both shows it's running.

@JeanLuc have you ever opened these files in a Windows machine? Could they have Windows line endings? Try running sed 's/\r//' patterns.txt > newFile and then tryingt he same approach with newFile. Also make sure that there are actually matches: remove the redirection, just run for chunk in pattern_chunk*; do echo $chunk;  grep -Ef "$chunk" target.fastq; done and see if there is any output.

Hello, I just did the sed newFile and re-run the while commands and it didn't work. Gave me the same result: echo is running but empty output.
I know for a fact there are matches, as by testing the first 10 lines of my patterns on this one specific target.fastq I get 500 hits. When I run grep normally like this, I get printed on screen all the hits highlighted etc (and I can save them to output easily, no issues there). But when I use the same patterns in a loop it gives me a blank output :/ grep just seems to not be doing it's job within the loops like when it's out of them

I would recommend not do this using grep. There is a program within the KAT suite called KAT filter for filtering reads by a set of kmers (that might require all your sequences to be of constant length, but I am not sure). Another software you might want to check out is Cookiecutter. These tools specialised in subsetting sequencing files are usually a lot more efficient than generic regular expression tools.

Hello! Thank you for your advice. I will try and see if the computer cluster possesses Cookiecutter already. Unfortunately, I don't necessarily have the required credentials to make any modifications (i.e. install and uninstall software) since it's an institution's cluster. But I will try and see nonetheless!

@JeanLuc you don't need administrative rights to install software, if you have a compiler (e.g. gcc) you can compile the source code and use it without moving the binary to shared directories (which is the part that requires administrative rights) or you can try to download the binaries and use them directly (https://github.com/ad3002/Cookiecutter/releases/tag/v1.0.0)


text:
Nextflow: How to initialize channel from CSV file but only after CSV file created?
Process A runs a python script which generates a CSV file in some absolute file path say "/mnt/x/blah.csv".
Process B needs to consume the CSV file line by line, after process A has finished.
The trouble is that the channel to input to B is initialized at the beginning before "/mnt/x/blah.csv" is created by A so it fails. This occurs even if I collect the result of process A as input to process B (because channel is initialized regardless)
nextflow.enable.dsl=2

def get_the_csv_records() {
    Channel
        .fromPath("/mnt/x/blah.csv")
        .splitCsv(header:true, quote: '\"')
        .map{ it.TheFieldIWant }
}

process process_a() {

    output:
    val("I'm done!")

    '''
    python.exe ./make_the_csv_file.py --path=/mnt/x/blah.csv
    '''
}

process process_b() {

    input:
    val(previous_step_is_done_dummy_variable)
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a()
    process_b(process_a.out.collect(), get_the_csv_records())
}

This fails immediately with "No such file: /mnt/x/blah.csv" because the channel is initialized before process_a has even started.
Same behavior occurs when using concat:
    process_b(process_a.out.collect().concat(get_the_csv_records()))

Question:

How can I make process_b() only read the input channel from the CSV file after it is created?

Alternatively, can I make process_b() trigger when process_a() is finished some other way?

Alternatively, can I add an intermediate step which reads process_a.out.collect(), then reads /mnt/x/blah.csv and outputs it into some existing channel which process_b() takes as input?

Another approach - what is the way to create a channel which outputs the results of splitCsv, when the file is not known at compile time?

Assuming you are using latest version of Nextflow which support the DSL 2 language, and assuming your python process output the file blah.csv, then the easiest way to move forward will be:
process process_a() {
    output:
    path("blah.csv")

    '''
    python.exe ./make_the_csv_file.py --path=blah.csv
    '''
}

process process_b() {

    input:
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a() \
        | splitCsv \
        | map{ it -> it.TheFieldIWant } \
        | process_b
}

You can also refer to the nextflow google group.

source: Stack Exchange

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

sort: votes


text:
Is there an efficient way to extract CIGAR strings for read pairs from bam files with python?
I am working with bam files and I have to check if reads of a specific position or their mates are soft clipped. So, I am looking for a fast way to extract the read pairs from a bam file in python. So far, I use pysam and fetch reads of a given position. From those reads, I have access to the CIGAR string and can check them for soft clipping. My problem is that I cannot access the CIGAR of the read mates. Is there a way to access information of read mates with pysam or another tool, especially CIGAR strings, without iterating over the whole .bam file?
Unlike in this question Extracting the CIGAR string from a BAM via Python? I am looking for the CIGAR of both reads of a pair.

I think that the easiest way is to work with query-name sorted groups of reads. In that case, mates will be adjacent in the sort and you can use that to extract the paired CIGARs.
If you are depending on position sorting to extract the alignments of interest efficiently, you can maybe do the following:

extract alignments in the region of interest plus a lot of padding to ensure that you are also getting the mates (this will obviously only work with standard short-insert shotgun libraries).

Query-name sort the list iterator using built-in sort() or sorted() with a key function that returns the query name.

Iterate over the pairs using adjacency in the iterator to extract paired CIGARs. Discard unpaired reads (first reads whose query name does not match supposed mates).

Code would look something like this (I have not run or tested this in any way!):
import pysam

def get_query_key(read):
  """read is pysam.AlignedSegment"""
  return read.query_name

def get_paired_cigars(bam, chr, start, end, padding=1000):
  """bam is pysam.AlignmentFile"""
  reads = list(bam.fetch(bam, chr, start-padding, end+padding))  # may fail at chr ends
  reads.sort(key=get_query_key)
  first_cigar = None
  mate_cigar = None
  cigars = list()
  qn = ""
  for read in reads:
    if first_cigar is None:
      first_cigar = read.cigarstring
      qn = read.query_name
    elif mate_cigar is None:
      if read.query_name != qn:
        first_cigar = read.cigarstring
        qn = read.query_name
      else:
        mate_cigar = read.cigarstring
    if first_cigar is not None and mate_cigar is not None:
      cigars.append((first_cigar, mate_cigar))
      first_cigar = None
      mate_cigar = None
  return cigars

I agree that I would like this to be simpler! Unfortunately BAM files are not terribly well-designed for using pairing information at that level of detail IMO, I have always had to use something like this. I would be interested to see a straightforward alternative to this approach from someone with more experience.

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


text:
How to activate an installed plugin (NetworkAnalyzer) on Cytoscape?
I need to use NetworkAnalyzer to analyse the nodes of a PPI network on Cytoscape, but the pre-installed app is not displayed on the Tools menu.
I checked the App Manager, but it shows "Installed" beside the app name and also it is not in the Currently installed apps list, which is why I'm unable to perform uninstall-reinstall to fix it.
I've once replaced the app by installing it manually, but despite of that there's no change in the Tools menu.
There's an option "Analyse Network" in the Tools, probably a substitute of NetworkAnalyzer, but it is not showing the usual attributes (Degree, BetweenessCentrality etc) which are displayed on using Network Analyzer app.


text:
How to create a DESeqDataSet and define experiment design before varianceStabilizingTransformation?
I have an RNAseq count matrix consisting of 2 groups (high, low) with 6 timepoints per group (T1,T2,...,T6) and 3 replicates per timepoint (rep1, rep2, rep3). So a 2-factor design with 36 samples in total (2x6x3=36).
I want to perform a varianceStabilizingTransformation (with blind=FALSE) of the matrix to later perform rhythm analysis (over the timepoints, but separately for each group). Now I wonder how to create the DESeqDataSet that the varianceStabilizingTransformation command needs as input and how to appropriately define the design formula associated with the DESeqDataSet.
Here a scheme of my count table (the numbers are random), which I would save in a .txt or .csv (named "matrix") and import with matrix <- read.table(FILE.txt, header=TRUE, row.names=1):
ID            high_T1_rep1   high_T1_rep2   high_T1_rep3   high_T2_rep1 ...
transcript1   1.23           1.45           1.67           1.89
transcript2   5.32           5.54           4.76           5.98
transcript3   3.22           3.44           3.66           3.88
transcript4   7.33           7.55           7.77           7.99
...

For the design I would now create a table (also txt. or .csv) as follows and import it with designdata <- read.table(FILE.txt, header=TRUE, row.names=1).
sample         group   timepoint
high_T1_rep1   high    T1
high_T1_rep2   high    T1
high_T1_rep3   high    T1
high_T2_rep1   high    T2
...            ...     ...
low_T5_rep3    low     T5
low_T6_rep1    low     T6
low_T6_rep2    low     T6
low_T6_rep3    low     T6

I went through several DESeq2 manuals, but still do not understand how to find the right definition for the design=... term when creating the DESeqDataSet. The varianceStabilizingTransformation seems to use the design information (if blind=FALSE), and thus I want to make sure, I provide it correctly.
From what I found out so far, my command for the DESeqDataSet would be something like this:
dds <- DESeqDataSetFromMatrix(countData=matrix, colData=designdata, design=???)

As options for the design I thought of e.g. "~ group + timepoint" or "~ group + timepoint + group:timepoint", but I have no idea how one or the other option would affect the following varianceStabilizingTransformation.
So while I have at least some rough idea what to do, my big concern is handling the data incorrectly and thereby distorting it without even noticing.
Any help you could give me is greatly appreciated. Thanks a lot in advance!

Are your real counts...counts?  Are they all whole numbers?  You realize that the vst values won't be used at all in DE calculations?  And what design you want depends on exactly what questions you want answered.

Thanks for your thoughts. These are TMM-counts (trimmed mean of M-values) I got from the salmon mapper and they have decimals. I also have raw counts (whole numbers), but those do not take differences in sequencing depth between samples into account, so I'm hesitant to use them.
I want use the vst values for a different analysis (R-package RAIN) and performing this transformation beforehand was recommended to me by a specialist on this method. In the analysis, rhythmicity is determined from the changes over the time points (T1, T2,...). This is done separately for the 2 groups (high, low).

The DESeq documentation is quite clear.  It wants raw counts.  Library size correction can be done on those in DESeq.

Thanks, then I will use the raw counts. Does DESeq automatically determine/consider sequencing depth differences between samples, or does info need to be provided? (at which step?)

Have you read the vignette?  Of course it will do a library size correction, or you can give it numbers to use, if you think your method is better.

I mostly used two pdf-manuals (Beginner's guide to using the DESeq2 package / Differential analysis of count data - the DESeq2 package), which seem similar to the webpage, which I guess you mean (Link). They do not mention the automatic correction for sequencing depth as clearly as the vignette, but now that I knew what to look for, I found it. I also saw that the correction is applied directly during the varianceStabilizingTransformation. It is however still not clear to me how to best define the design term of the DESeqDataSet so the later rhythm analysis works ok.


text:
How to run a phylogenetic analysis with non-independent data?
These are the facts: I am working with some morphological features. Each trait can have one of the following three states: 0, 1, 2.
State 0 is always independent.
But if there is a trait with state 1, then there must be, in another column of the same row, a state 2 (and vice versa, obviously).
The traits can't be encoded separately because, eg, in the row_01, at the trait_01, the state 1 is ligated to the trait_20 with state 2. But in another row (row_02), the trait_01 with state 1 could be ligated to another trait (trait_33) with state 2. My data behaves like this because I'm working with morphology in... molecules.
So far, I know that the data of an array in phylogeny must be independent, so I am wondering if any of you know of any way to approach this problem (paper, algorithm, your own propose, etc.).
solving conditions:

I can't discard the traits with state 1 or 2.

Cheers!

I am not sure that I understand. If the traits fully determine each other, why are they encoded separately? It might help to see a sample of the data. If traits were even slightly decomposable one might suggest an ordination method (like PCA) to find independent eigentraits, for example http://adn.biol.umontreal.ca/~numericalecology/Reprints/Guenard_et_al_MEE_2013.pdf.

Since I last looked at this question, the editing appears to have obfuscated it. e.g. My data behaves like this because I'm working with morphology in... molecules. what does it mean for data to be "ligated"?


text:
How to extract subset of protein structure (PDB format) file based on a subsequence of that protein
I looking at a particular protein structure called 2LY4 accessible from RSCB PDB website.
The corresponding fasta sequence for that structure is this:
>2LY4_1|Chain A|High mobility group protein B1|Homo sapiens (9606)
GKGDPKKPRGKMSSYAFFVQTCREEHKKKHPDASVNFSEFSKKCSERWKTMSAKEKGKFEDMAKADKARYEREMKTYIPPKGE
>2LY4_2|Chain B|Cellular tumor antigen p53|Homo sapiens (9606)
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPL

And the PDB format file can be downloaded here (1.7MB). The entire is to large to be pasted here.
What I want to do is to extract the subset of PDB format based on the subsequence in fasta above.
Namely Chain A starting  from 1st residue to 30th residue.
GKGDPKKPRGKMSSYAFFVQTCREEHKKKH

How can I do that in R or Python?

In general, if you simply want to extract that part of the PDB file, you could loop over it (it's plain text) and check the fields you're interested in:
with open('2ly4.pdb') as pdb:
    for line in pdb:
        if line[:4] == 'ATOM':
            chain = line[21]
            res_idx = int(line[23:26])
            if chain == 'A' and 1 <= res_idx <= 30:
                print(line[:-1])
        else:
            print(line[:-1]) # alternatively, skip with continue

However, this doesn't quite work with the file you've linked because it actually contains ten different conformations (models). This is better handled with a proper parser like prody:
import prody

model = prody.parsePDB('2ly4.pdb', model=1)
extract = model.select('chain A resnum 1to30')
prody.writePDB('extract.pdb', extract)

This will also allow you to do more complex atom selection.

Neat! The downside is that it removes all the headers including important ones, such as CRYST1.

If you remove model=1 all models are written
(but with a small change: in MODEL 3 the charge -1 that in the original file is on atom 390 (OE1 in Glu 25) is moved to atom 391 OE2).

Right, parsing the file by hand gives you more control over the output, at the expense of reinventing the wheel.

Here is an example how to do it using another Python library with another selection syntax.
Requires pip install gemmi.
import gemmi

st = gemmi.read_structure('2ly4.pdb')
subset = gemmi.parse_cid('/*/A/1-30').copy_structure_selection(st)
subset.write_minimal_pdb('output.pdb')

As @marcin suggested, it's probably better not to rely on ad hoc parser and use functionalities already provided by specialized libraries.

If you need to do it only once, for this one PDB file,
it's easier to do it manually.
When you open the file you can see it has 10 models from NMR.
Do you want to extract a subset of one model or of all 10?
In this file the residues are numbered according to the sequence.
So you find the line where residue 31 begins in the first model:
ATOM    500  N   PRO A  31      18.844  -9.699  -4.620  1.00  0.00           N

and delete all lines until
ENDMDL

Then you do the same for other models, or you just remove other models, depending on what you need.
You may also want to remove headers – from the beginning of the file until the CRYST1 record.


text:
How do I increase the sensitivity of Salmon Alevin?
I'm in a little bit of a bind with targeted single-cell sequencing. I'm trying to match up our reads to the targeted amplicon panel (418 targets), and all but one have matched successfully with Salmon Alevin. I have a transcript that produces many 75bp reads that only match the first 28bp of a target, and Alevin is ignoring those reads as a match (the 28bp sequence is genome-unique).
Here's an example read, with matched sequence:
qual FFFFFFFFFFFFFFFFFFFFFFFFFFFF,,,::,,:,:,,,,,,,,,,,,,,,,,:,,,,,,,,:,,,,,:,:
read GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACGCGCGCCGCGCCCCCCCGCCCCCGGGGGGGGGGGGGGG
SNPs ...................................G.G.G..G.GC.....CG...CCGGGGGGGGG..GG..
 ref GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACCCCCACCCCCACCCCCACCCCATTCTTAAATTGTTTGG

As you may be able to see from the above sequence, the reads for this particular transcript get up to 28bp (sometimes with one or two errors), then the wheels start coming off the sequencer, showing an abundance of Gs. This was a NovaSeq run, and a polyG read is the same as no signal using its two-dye process, so I suspect that the sequencing stops after ~35bp, rather than actually synthesising Gs from thin air.
I have lots of these reads; there are about 400k reads that match the reference transcript (Cxcr4, just in case it matters), and about 250k of those reads have a G homopolymer that's at least 8bp in length.
Any ideas on what to do?
FWIW, the above match was taken from a Bowtie2 alignment I did on the reads. I could use Alevin for cell barcode correction only and do the UMI deduplication manually, but it seems a shame to do that when it's only one transcript (from 418) that's spoiling a fully-automated process.
I did try adding in an additional transcript that included the full sequence of a short-match read, but there seem to be too many errors in the non-matching sequence for any reads to count (at least... I assume that's why the counts are not coming through; I'm not quite sure how I can work out what's going on behind the scenes).
Is there any way to tweak Salmon Alevin to get it to accept these shorter matches? If I can change the kmer size, how do I do that?
End trimming is unlikely to help, because that's going to shorten the reads and make them less likely to map to the target sequence. I want a way to include these reads as valid hits to the target sequence.

So, I think there are a few potential options here.  Alevin is using selective-alignment internally to determine the mappings for the reads.  So, even if you have a k-mer supporting the mapping, if you have gibberish for the rest of the read, the mapping score is going to be very poor and you're not going to recover that mapping locus.  There are some parameters you can mess with to try to increase the sensitivity of mapping, though you should always be aware that upping the sensitivity a lot will mean there is less aggressive filtering early in the mapping process which can slow things down.  The main parameter you should tweak to allow lower alignment scores is --minScoreFraction; the smaller this value the lower (worse) will be the permitted scores for valid mappings.  Additionally, you can combine changes to the --minScoreFraction with the --softclip option to simply score the non-matching end of read as 0 (rather than a long deletion or large set of mismatches); this will lower the effect of bad read ends.
Another alternative you might consider is to try out alevin-fry for processing this data, and mapping the reads with the --sketch flag.  This will only check certain structural constraints of the mapping and won't force a particular minimum alignment score.  The caveat here is that dropping the requirement for alignment scoring can increase the rate of spurious mappings (reads deriving from elsewhere that spuriously map to some indexed target simply by virtue of sharing some sequence).  In transcriptome wide analyses, this can be greatly mitigated by building an appropriate index and using a UMI resolution strategy that complements this expanded index. We've not examined this in detail for targeted single-cell sequencing.  However, if you want to give it a try and run into any issues, feel free to drop an issue over on GitHub and we'd be happy to help work through it.

Thank you; tweaking the minScoreFraction did the trick. I went from a max Cxcr4 read count per cell of 1 with default parameters, to 4 with minScoreFraction=0.5, to 26 with minScoreFraction=0.33.

@ATPoint commented that the default kmer size in Salmon indexing is 31, which led me to a solution that involved altering the index generation step:
Instead of:
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_${targetName}

I now have:
salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_k11_${targetName} -k 11

In other words, I changed the kmer size from its default of 31 to 11, which matches the BD-recommended seed substring length for the standard Bowtie2 approach ("-L 10", but I increased that to 11 because Salmon needs odd kmer sizes). I'm pretty sure that my 28bp-mapping reads won't be picked up with a kmer size of 31. Because this is targeted sequencing, I don't expect that a reduced kmer size will be an issue, but I'll keep a watch out for that.
I'm currently re-running the analysis. If it still doesn't work, then I'll see about adding in the tweaks suggested by @nomad.
Update: still not working, so I'll shift onto score tweaks.


text:
How to extract docked ligands from .maegz docking results file sorted by glide gscores in Schrodinger Maestro?
My VLS run completed yielding millions of ligands and I would like to extract top-scored 20000 ligands in a separate file. Incorporating the .maegz output in Maestro is not an option as it fills up all the memory and keeps crashing. Are there any command-line options for this task?
PS: I tried using rdkit in python to parse the output file after converting it to .sdf, however, rdkit skipped some ligands that it could not read.

It sounds like you are saying your file is already sorted? Is this correct?
RDKit Mae supplier is an iterator as is the gzip reader, so does not load all in memory, so that is the best option given that loading all 1e6 into memory crashes your system. What is the fraction of unreadable compounds in RDKit? And does it matter? These are likely failing sanitisation —the sanitize=False flag. If for some reason none can be ignored, circumventing this is hard. The removeHs=True flag is rarely/never the cause even though the description of several aromatics require protons, e.g. an indole ring is c1ccc2c(c1)cc[nH]2, the proton removal happens after loading. If you have an identifier and a list of SMILES you can correct dodgy bond orders etc. But making an unsanitisable molecule sanitary is tricky and rarely worth it as the output molecule is likely rubbish anyway.
Plus loading in RDKit will allow you do some additional calculations, such as ligand efficiency (best for sorting anyway), partition coefficient etc etc. As opposed to simply slicing your sorted file.
Alternatively, I believe Open Babel does not load all into memory and has a splitinto argument to make the file smaller and can convert gzipped Mae file to SDF.


text:
Total pairs of amino acid substitution
I am reading a blog and it says:

The numbers for identities and replacements used for calculating the overall alignment score in the expression above are usually presented in the form of a 20 x 20 matrix (20 is the number of the most common amino acids). In total there are 210 possible replacement pairs (residues replacing each other) of amino acids, which includes 190 pairs of different amino acid substitutions + 20 pairs of identical substitutions

My question is how are there 210 replacement pairs. What is the formula being used here? Insights will be appreciated. The link to blog is :
https://proteinstructures.com/sequence/amino-acid-substitutions/

There are 190 pairs of different amino acids:
$$ {20\choose 2} = 190$$
and 20 pairs of identical amino acids (AA, TT, WW, etc). Combining these, you get all 210 pairs referenced in the article. Another way of seeing this is dividing the 20x20 matrix into 3 parts: the diagonal, all values above the diagonal, and all values below. We have 20 values on the diagonal with 380 values to split evenly across the upper and lower parts, leaving 190 for each. Therefore, the lower (or upper) part plus the diagonal comes out to 210.

Thank you. I understand it now.


text:
Network analysis: confidence score (STRING database) validation
I am working on network analysis of differentially expressed genes, but the problem I am facing is that the gene set is very small (around 16 genes only), so when I'm trying to build an interaction network with these genes using STRING, the nodes starts to dissociate from the interaction when I set the confidence score to 0.400 or above (as the genes are related to Homo sapiens).
DE gene network has to be made in cytoscape, for that STRING ppi interactions are needed, that interaction network will tell us co-expression, homology etc and provide us a combined score.
So, to build an interaction network between these genes, I have to decrease the confidence score to 0.150 since it is only with this score that I am finding interactors.
So, I want to ask whether this low confidence score is valid for this short gene set related to Homo sapiens or, if not, what should I do in this case?

Why do you think there should be a connection between these genes or proteins? Are you working with genes and trying to build a gene regulatory network or are you working with proteins and trying to build a protein-protein interaction network?

We really can't help unless you explain why you think there should be connections between these genes and why that connection would be detectable from PPI networks as opposed to gene regulatory networks.


text:
Programmatically retrieve Metadata from SRA Run Selector
I previously asked a question about how to retrieve the Accession List associated with a SRA project. The answer was:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format runinfo
where PRJNA491191 is the bioproject that I'm interested in in this example. How can I do the same for the metadata file (the arrow in the picture below is pointing to accession list. It should point to metadata, which is what I'm interested in now)?

Depending on exactly which information you seek, you can try:
esearch -db sra -query 'PRJNA491191[bioproject]' | efetch -format native -mode xml

This will return an XML with additional information that is not present in the runinfo table. You can use the Entrez Direct tool xtract to parse the output XML and extract the fields of interest into a tabular format.
Alternatively, you can try:
esearch -db sra -query 'PRJNA491191[bioproject]' | esummary

to retrieve a summary output that includes fewer fields than the efetch option shown above. Perhaps all of the fields that are of interest to you are present in the summary output... it will save you time and effort in parsing the XML and retrieving the information you want.
For example, you can use the following to retrieve some fields using xtract from the esummary output:
$ esearch -db sra -query 'PRJNA491191[bioproject]' \
  | esummary \
  | xtract -pattern DocumentSummary -element Summary/Title,Statistics@total_size,Study@acc
GSM3391430: hsf-1(OE) rep3; Caenorhabditis elegans; RNA-Seq  731402855  SRP161795
GSM3391429: hsf-1(OE) rep2; Caenorhabditis elegans; RNA-Seq  703822418  SRP161795
GSM3391428: hsf-1(OE) rep1; Caenorhabditis elegans; RNA-Seq  715482490  SRP161795
...

Thank you. So I want to create a csv file from the output. Also, I'd like to keep all columns. How can I do this? I'm trying esearch -db sra -query 'PRJNA491191[bioproject]' | esummary | xtract -pattern DocumentSummary -sep "," -element but I'm not sure what should go after -element since I want to keep everything

Also, I can't specify all columns because I don't know what they are beforehand (I know I can look them up from the output but I'm doing this for many projects so I practically can't do this manually)

You have to specify the elements to be able to retrieve the fields you are interested in. Unfortunately, you cannot retrieve everything using xtract without knowing the field names. I will update my answer above with an example.


text:
How to convert CRAM file with 10x data in three fastq files
10x Genomics data are stored in three FASTQ files, besides the standard R1 and R2 reads, there is also a I1 file with some metadata. Sometimes however they are shipped in a single bam/cram file (e.g. the data from the Darwin Tree of Life). How can I convert the BAM or CRAM files to FASTQ?
-- edit --
Apparently, the I3 file is not 10x specific thing. It is an index generated when demultiplexing any Illumina data. See this question for a nice explanation.

Isn't I1 the index read not unique to 10x data? But normally demultiplexing Illumina data with bcl2fastq does not generate these reads with default arguments

@Chris_Rands (after a bit of googling) you are right. I just not aware of its existence before I run into this question on the tDToL GitHub. And then I thought I might as well put it here so it's easier to google it.

All can be done with samtools. This is how the Darwin Tree of Life folks convert it:
samtools fastq -@4 -i \
  -1 ${sample}_S${tag}_L%03s_R1_001.fastq.gz \
  -2 ${sample}_S${tag}_L%03s_R2_001.fastq.gz \
  --i1 ${sample}_S${tag}_L%03s_I1_001.fastq.gz \
  --index-format i8 ${lane}.cram

I don't know if that solution will work with 10xGenomics files.  You have to recreate R2 out of the contents of tags; there will be no reads with the read2 binary flag set.

@swbarnes2 I just have seen this thread: https://github.com/darwintreeoflife/darwintreeoflife.data/issues/2#event-4859159249 and tried to make it more googlable (because I do think there will be quite a few people who will be downloading the tDToL data).

I saw that thread too.  I don't think it's right.  samtools fastq -2 will work fine for bams from many applications, but not 10XGenomics.

@kamil, can you post a couple of lines of the BAM file showing how tags for R1 and R2 appear in 10x data? I can probably quite easily add support for this in https://genozip.com/sam2fq.html

You'll want to use the 10x provided bamtofastq tool to preserve the indices properly: https://github.com/10XGenomics/bamtofastq


text:
PCA plot in R coloured by sample type
I'm a biologist, not a programmer so please be gentle.
So I have a dataset that looks like
Genes  Patient1   Patient2   Patient3
A          324      433         343
B          431       342        124
Z          232       234        267

then I have the sample sheet where it contains sample info like:

Patient1 - Healthy
Patient2 - Disease
Patient3 - Healthy
(110 samples)

I am using:
library(ggfortify)
df <- dataset
pca_res <- prcomp(df, scale. = TRUE)

autoplot(pca_res)

Then I want to do
autoplot(pca_res, data = ?, colour = '?')

I wish to use the info from the sample sheet to colour my PCA based on the state (healthy/disease) using the autoplot function. Is there a way to do this?

I think you need to provide more information for a minimal reproducible example.

Given that you probably do not want to invest lots of time in such a simple task I'd just read the PCAtools vignette over at Bioconductor. That package can probably do all you want including making pretty plots, and the heavy work is coded under the hood so for the end user it comes down to running a couple of commands in R.

I've found this the way to work with autoplot (ggfortify's prcomp autoplot) - Use a df that has the field you need for color as the first column (to make things easy), then do:
autoplot(prcomp(df[,-1], data=df, color='column1_name')

If the color-by field is not the first column, do:
autoplot(prcomp(df[,!(colnames(df) == "color_by_column_name")], data = df, color = "color_by_column_name")

EDIT:
That being said, you're going to need to process your data so patients are along rows instead of column. Transpose your dataset so rows are columns and columns become rows; change patient ID to rownames and Gene name to colnames, then merge the patient category dataset with the expression dataset. Once you merge, your dataset should look like:
         Category    A      B      Z
Patient1 Healthy     324    431    232
Patient2 Disease     433    342    124
Patient3 Healthy     343    124    267

And then,
autoplot(prcomp(df[,!(colnames(df) == "Category")]), data = df, color = "Category")

Hi Ram,Thanks a lot for taking the time to answer me. I did exactly what you told me above and there is an error
Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric
Any ideas why is that ? I am not sure at all

Ensure that the columns you're passing to prcomp are all numeric. That's probably what's going wrong here.

sorted, thanks a lot :)

Can't say I've ever used autoplot, but this is fairly easy to achieve using "base" ggplot.
One thing autoplot does is put the axis on the same scale reflecting the percentage variation captured by each PC - I think this is a good thing, but most PCA plots you'll see don't do this.
As you haven't given us any info on your metadata, I'll just make a few assumptions - lets make it a data.frame called "colData" with a "Status" column.
prcomp stores the PC scores in "pca_res$x"
# generate some metadata
colData <- data.frame(row.names=rownames(pca_res$x),status=as.factor(sample(c("diseased","healthy"),nrow(pca_res$x),replace = T)))

# get % variation for each PC
pca_res$percentVar <- pca_res$sdev^2/sum(pca_res$sdev^2)

# change scale to reflect % variation
d <- t(data.frame(t(pca_res$x)*pca_res$percentVar))

# merge metadata and modified PC scores by row.names (by=0)
d <- merge(d,colData,by=0)

# make the base plot (PC1 vs PC2)
g <- ggplot(d,aes(x=PC1,y=PC2,colour=status))
# add the points
g <- g + geom_point(size=2)
# add a "wow" factor to the points
geom_point(size=4,alpha=0.5)
# display the plot
g

There may be a way to do this with autoplot - e.g.
autoplot(pca_res,colour=as.numeric(colData$status))

gets you part of the way there, but I think using the standard ggplot syntax is easy enough anyway.


text:
How to identify the DNA sequence of a gene in the complete genome sequence FASTA format file?
I need to identify the sequence of a gene in the complete genome sequence . I thought it was simple, instead it is not a straightforward task !
My method was the following: I downloaded the FASTA format file of the complete genome sequence of the archaea organism that I am studying, that is CP003685, and then I downloaded the CDS file from NCBI by clicking on ">Send to>Coding Sequencing>FASTA Nucleotide>Create File". In this way I have a txt file (always FASTA format file) in which I find a sequence for each gene (if I am not wrong).
So for example, I should find, in the complete genome sequence, the sequences of the CDS file (the last that I downloaded), right?
For example, if I consider:
>lcl|CP003685.1_cds_AFN03983.1_1007 [locus_tag=PFC_05180] [protein=DNA protection protein DPS] [protein_id=AFN03983.1] [location=complement(940023..940580)] [gbkey=CDS]
ATGCCAGAGCATAATAGGAGATTAGTTGAAAGGACAGGAATAGATGTAGAGAAGCTGTTGGAACTTCTCA
TTAAAGCGGCTGCAGCAGAGTTCACAACGTATTATTACTACACTATACTTAGAAACCATGCAACAGGTCT
GGAAGGAGAAGCAATTAAGGAGATTATCGAAGATGCTAGACTTGAGGACAGAAACCATTTTGAAGCCCTC
GTGCCAAGAATATACGAGCTTGGAGGAGAATTACCAAGGGATATTCGGGAATTTGCAGATCTAGCTTCTT
GTAGAGATGCCTACCTACCAGAAGAGCCCACTATAGAGAACATTCTAAAAGTCCTACTCGAAGCTGAGAG
ATGTGCAGTTGGAGTCTATACAGAAATATGCAACTACACTTTTGGAAAGGATCCAAGGACTTATGATTTG
GCCCTTGCCATATTACATGAAGAGATAGAGCATGAAGCATGGTTTGAAGAGCTCTTGACTGGAAAGCCGA
GTGGTCACTTTAGAAGAGGAAAGCCTGGGGAAAGTCCATATGTCTCAAAGTTTTTGAAAACTAGATAG

I should find in the complete genome sequence, this sequence, right ? But this is not. Why?
I suspect that I miss very important passages.
Thank you in advance.

See How to identify genes from a genome assembly of C. Elegans?. Exonerate and wise2 should be able to do this.

I simply copy the sequence of the gene (like the one I have written in the post) and in the file with the complete genome sequence (CP003685) I click ctrl+F and paste the reverse complement of the sequence of the gene and search it. But it does not find it. Sorry if I did not really specified this very crude way to do it :') . For the introns I have no idea...

It is a non-professional way but I do not understand why it does not work since the two sequences should be the same... the two files I am considering (as said below) are only two different ways to represent the genome DNA sequence, right ?

Well yes, kinda, but you have line breaks. You will never find the entire thing unless you convert it to a one-line format first. Try searching for a smaller section, say ~10 nucleotides and if the first one doesn't work try some others. You need to find a sequence that appears in a single line otherwise you won't find a match.

Thank you very much @terdon. I did not know this and now it seems to be finally clear putting all pieces together !

I would just blast it. When blasting locally, you need to first make a database from your genome, so assuming you got the command-line version of blast installed you can do something like
makeblastdb -in my_study_genome.fa -dbtype nucl
blastn -max_target_seqs 10 -db my_study_genome.fa -evalue 1e-10 -outfmt 6 -query my_downloaded_gene_of_interest.fasta -out gene_in_my_study_genome.blast

blast search for anything that is "similar enough" and also considers both strands, so for example, if there is an SNP in the middle of the gene. Hence even if the searched gene and the reference genome are from different sequencing runs / individuals, blast will still find it. Blast will also handle linebreaks for you (in fasta files the nucleotide sequence is quite often wrapped by 60 or 100 characters a line). Long story short, there is not much that can go wrong with blast (while with full text search you can easily shoot yourself in foot).
-- edit 2 --
As @terdon suggested, potentially if you would locating a transcript sequence (excluding introns), you would do better with software like exonerate or genewise, because they model exon/intron boundaries. However, even in those cases blast will give you an approximate position of the gene.

@Manuela this really should be done using a tool that models genes (try exonerate or genewise). Blast is great but will not be able to model introns correctly if your species has splicing.

In your sequence header you will see:
[location=complement(940023..940580)]

It's a reverse strand gene. Reverse complement the sequence, then search the genome sequence.

Thank you very much for the answer but it is still not working . I do not understand why. I used this software to make the reverse complement, copied and searched in the complete genome sequence but I still do not find it.


text:
Find selective markers for one cluster automatically
I detected all markers for each cluster and checked 20 top markers for my clusters one by one. Is there any way that I figure out the selective markers automatically from all the clusters?
For instance, in seurat tutorial, I found comparing one cluster two or three other clusters, but I want to see for instance what are the 10 top markers or genes that are expressed only in cluster of interest.

Directly copy-pasting from one of the Seurat vignettes:
# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
    group_by(cluster) %>%
    top_n(n = 2, wt = avg_log2FC)

Change n to 10 in the top_n() call and filter for your cluster of interest and you will get what you want.

Thanks a lot,  I realized this code exist in seurat web page, but I didn't know that does the same. Thanks for answering.


text:
Renaming samples in vcf file
I am currently asking myself how to rename samples in a vcf file. This is necessary, because the naming of several samples included in thr vcf file was unfortunately not be chosen by me.
I am loading the vcf file via readVcf() from the VariantAnnotation package in R.
To check the sample names, I can use samples(header(vcf.file)), which will give me the names as characters.
It looks like this:

and so on..
(I am actually very sorry that I cannot give a reproducible example here :( )
I also have a vector with all the new names. For example
new.names <- c("1","2","3","4")

I tried to exchange the "old names" with the new.names vector:
samples(header(vcf.file)) <- new.names

which leads to the error
Error in samples(header(data.vcf)) <- names.sample :
  could not find function "samples<-"

I also tried to use the rutilstimflutre package, which has a already designed function for this, but unfortunately the package is outdated and not available for the newest R version.
Is there a way to do this? Perhaps a better way to assign the new names to the old ones?
Thanks a lot.

Easiest to use bcftools, if samples.txt has the new sample names:
bcftools reheader -s samples.txt -o new.bcf old.bcf

Thank you, since i only have a brief knowledge of R, is there an R solution?

Maybe - I'm not sure. This approach doesn't require any knowledge of any programming - just doing this should work perfectly. Generally I prefer not to use R for working with .vcfs because it's slow to load them into memory. If you are working with .vcfs I would highly recommend getting familiar with bcftools as it is very fast and memory efficient.

This proved to be so simple that it took a while to discover: Just use colnames.
x = readVcf("in.vcf")
colnames(x) <- c("new1", "new2")
# At this point, colnames(x) and samples(header(x)) will show different names

# ... but if you write it out, it should use the colnames for both:
writeVcf(x, "out.vcf")
y = readVcf("out.vcf")

# samples(header(y)) and colnames(y) should both be "new1","new2" here

Thank you for posting a solution in R.

Here is an example in python how to rename samples using pyvcf submodule I wrote:
>>> from fuc import pyvcf
>>> data = {
...     'CHROM': ['chr1', 'chr2'],
...     'POS': [100, 101],
...     'ID': ['.', '.'],
...     'REF': ['G', 'T'],
...     'ALT': ['A', 'C'],
...     'QUAL': ['.', '.'],
...     'FILTER': ['.', '.'],
...     'INFO': ['.', '.'],
...     'FORMAT': ['GT', 'GT'],
...     'Steven': ['0/1', '1/1'],
...     'Sara': ['0/0', '0/1'],
... }
>>> vf = pyvcf.VcfFrame.from_dict([], data)
>>> # vf = pyvcf.VcfFrame.from_file('your_file.vcf')
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT Steven Sara
0  chr1  100  .   G   A    .      .    .     GT    0/1  0/0
1  chr2  101  .   T   C    .      .    .     GT    1/1  0/1
>>> vf.df = vf.df.rename(columns={'Steven': 'John', 'Sara': 'Rachel'})
>>> vf.df
  CHROM  POS ID REF ALT QUAL FILTER INFO FORMAT John Rachel
0  chr1  100  .   G   A    .      .    .     GT  0/1    0/0
1  chr2  101  .   T   C    .      .    .     GT  1/1    0/1
>>> # vf.to_file('renamed.vcf')


text:
Positive marker and negative marker in scRNA-seq datasets
I am learning about the marker genes and clustering by Seurat in scRNA-seq datasets. However, I am confused about the term positive marker and negative marker. What I understand is the marker genes of a cluster is the genes which are highly expressed (high counts) in the cells of that cluster. Let's say cluster A has 200 genes. My questions are:
Can I say that those 200 genes are marker genes for cluster A?
How to determine which marker genes are positive marker genes and which are negative marker genes?

In scRNA-seq clusters are groups of cells with similar transcriptomes typically defined by graph based methods (e.g. leiden, louvain) or more rarely k-means.
Marker genes could be genes that are known to be highly expressed (or better specifically expressed) in a particular cell type (e.g. CD3E is a T cell marker). You can visualize the expression of such marker genes on a UMAP/t-SNE without doing any clustering.
When talking about positive and negative markers, I assume the reference is too differentially expressed genes between different clusters. Typically simple differential expression tests are conducted (e.g. Wilcoxon rank sum) between each individual cluster and the other collective clusters. The upregulated differentially expressed genes could be considered positive markers and downregulated genes negative markers.

What do you mean upregulated and and downregulated? I do not understand about them? Also, the 200 genes I mentioned, basically we cannot say those genes are marker genes for cluster A?

A positive marker is a gene that is overexpressed in a cell population compared to other populations.
A negative marker is a gene that is underexpressed in a cell population compared to other populations.
I think the Seurat FindMarkers function has an option to output Log2FoldChange as well as the marker names. In this case the positive markers will have a positive Log2FoldChange, and the negative markers will have a negative Log2FoldChange.
Both positive and negative markers can be used to describe a cell population. In the research that I've done in the past, there's a particular type of dendritic cell that is only distinguished from another dendritic cell type based on cell surface markers that it doesn't have; our own name for these DCs was "Triple negative".


text:
Extracting WDL map keys as a task
I want to get the array of keys from a map data structure, so I can perform array-based functions against the map in subsequent tasks. I am limited to draft-2 Workflow Description Language (WDL), so I don't have access to the newer spec's extraction functions. I thought I could achieve this by writing a task:
task extract_keys {
  Map[String, String] map

  command {
    cut -f1 ${write_map(map)}
  }

  output {
    Array[String] keys = read_lines(stdout())
  }
}

This task works just fine and, as far as I can tell, returns the array of map keys. However, when I try to feed this array into subsequent tasks, I get an error:
workflow test {
  Map[String, String] my_map = {
    "a": "foo",
    "b": "bar",
    "c": "quux"
  }

  call extract_keys { input: map = my_map }

  scatter (k in extract_keys.keys) {
    call task_foo { input: k = k, s = my_map[k] }
  }
}

Call input and runtime attributes evaluation failed for task_foo

I don't have access to the full error message or the logs, unfortunately. The frontend I'm using hides these away.
(Incidentally, I'm aware that you can scatter over maps, although it's undocumented, per https://github.com/openwdl/wdl/issues/106#issuecomment-356047538. This code is just illustrative; I wish to do other things, such as length and cross, that only accept arrays.)
Is it not possible to scatter over the output of a previous task?

I seem to have solved this myself by virtue of the linked undocumented feature (scattering over a map). It turns out you don't need a separate task at all. This works:
workflow test {
  Map[String, String] my_map = {
    "a": "foo",
    "b": "bar",
    "c": "quux"
  }

  scatter (pairs in my_map) {
    String keys = pairs.left
    # String values = pairs.right
  }

  scatter (k in keys) {
    call task_foo { input: k = k, s = my_map[k] }
  }
}


text:
How do I remove allele annotations from SNP Ids in .bim file?
Hello every one i need help in editting my .bim file. So my .bim file looked like following. I want to convert the ids with chromosomal position from chr1_867635_C_T to chr1_867635 while others remained as it is. I want to remove all the alleles from my SNP IDs, not just specific ones.
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

The simplest approach is to remove the first occurrence of a _ followed by a capital letter, another _ and another capital letter from each line:
$ sed 's/_[A-Z]_[A-Z]//' file
1       rs2880024             866893  A       G
1       chr1_867635       867635  A       G
1       chr1_869303       869303  A       G
1       chr1_873558       873558  C       A
1       exm2275405           879439  A       G

Alternatively, find the first occurrence of chr followed by numbers, then an underscore and more numbers and remove all non-whitespace characters after the last numbers (note that this assumes GNU sed, the default on Linux systems):
$ sed -E 's/(\schr[0-9]+_[0-9]+)_\S+/\1/' file
1       rs2880024             866893  A       G
1       chr1_867635       867635  A       G
1       chr1_869303       869303  A       G
1       chr1_873558       873558  C       A
1       exm2275405           879439  A       G

If your file is more varied than you show, you can use this to ensure you only make the change on the 2nd field:
$ awk -F'[ ]' '{ sub(/_[A-Z]_[A-Z]/,"",$2)}1' file
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

Or, the same but only matching A, C, T or G:
$ awk -F'[ ]' '{ sub(/_[ACTG]_[ACTG]/,"",$2)}1' file
1       rs2880024             866893  A       G
1       chr1_867635_C_T       867635  A       G
1       chr1_869303_C_T       869303  A       G
1       chr1_873558_G_T       873558  C       A
1       exm2275405           879439  A       G

The output of these commands should be piped out to another file (e.g. using >)

Yep. Either that, or add an -i to the sed commands and a -i inplace to the awk (assuming GNU awk) commands to edit the file in place.


text:
What is the correct way of setting the genotype after splitting multi-allelic sites in a VCF file?
This question was also asked on Biostars
Consider this minimal VCF file which contains a multi-allelic call made by a somatic variant caller on a human tumor sample:
##fileformat=VCFv4.3
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##FILTER=<ID=clustered_events,Description="Clustered events observed in the tumor">
##FILTER=<ID=multiallelic,Description="Site filtered because too many alt alleles pass tumor LOD">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  foo
chr10   43613843        .       GC      TC,CC,TT        .       clustered_events;multiallelic   DP=5180 GT:AD:AF:DP:F1R2:F2R1:SB        0/1/2/3:0,5154,11,15:0.994,0.00229,0.00306:5180:0,2507,4,8:0,2576,6,5:0,0,2585,2595

So, we have a site with three ALT alleles and a genotype of 0/1/2/3.
The VCF specification has this to say about the genotype field:

GT (String):  Genotype, encoded as allele values separated by either of/or|.  The allele values are 0 for the reference allele (what is in the REF field),  1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on.

As I understand the VCF specification, given the ALT field TC,CC,TT, the genotype of 0/1/2/3 is telling me that the reference sequence is absent but there are heterozygous calls for each of the three variants in ALT. This makes sense because this is a tumor sample, so presumably what we have here is three populations of cells, some with the GC=>TC change, some with GC=>CC and some with GC=>TT. Given the very high support for the GC=>TC change, I would assume this is actually a germline variant and there are small populations of tumor cells in my sample with either the GC=>CC or the GC=>TT variant. Alternatively, and quite likely, the latter two variants are sequencing artefacts.
Regardless, I need to split these into one variant per line. As I understand the specification, after splitting, I should be left with three variants and these specific genotypes (I am showing only the genotype here, for simplicity):
chr10   43613843   .   GC    TC   GT 1/0
chr10   43613843   .   GC    CC   GT 1/0
chr10   43613843   .   GC    TT   GT 1/0

Is this actually correct? Unfortunately, standard tools like bcftools produce something different:
$ bcftools norm -m -any foo.vcf 2>/dev/null | grep chr10 | awk '{print $1,$2,$3,$4,$5,$9,$10}'
chr10 43613843 . GC TC GT:AD:AF:DP:F1R2:F2R1:SB 0/1/0/0:0,5154:0.994:5180:0,2507:0,2576:0,0,2585,2595
chr10 43613843 . GC CC GT:AD:AF:DP:F1R2:F2R1:SB 0/0/1/0:0,11:0.00229:5180:0,4:0,6:0,0,2585,2595
chr10 43613843 . GC TT GT:AD:AF:DP:F1R2:F2R1:SB 0/0/0/1:0,15:0.00306:5180:0,8:0,5:0,0,2585,2595

Here, we have sites with a single ALT variant but which have multiple genotypes and that doesn't make sense to me. How can we have a 0/0/0/1 genotype and only one ALT allele? Is bcftools assuming this is a triploid genome?  Is this a bug in bcftools or am I misunderstanding the VCF specification? What is the correct way to display the genotype after splitting this kind of multi-allelic call?


text:
Error: "fewer tokens than expected" while using --update-name flag in plink
I am trying to update my rsIDs in .bim file by using --update-name flag in plink but getting this error message. what could be the reason of it?
plink --bfile updated_ff  --update-name updated_ref.csv --make-bed --out test_updatedrs

PLINK v1.90b6.21 64-bit (19 Oct 2020)
Options in effect:
  --bfile updated_ff
  --make-bed
  --out test_updatedrs
  --update-name updated_ref.csv

Hostname: login1
Working directory: /gpfs/home/bsc08/bsc08856/phs001273_dbgap
Start time: Thu Jun 10 14:23:44 2021

Random number seed: 1623327824
96310 MB RAM detected; reserving 48155 MB for main workspace.
Allocated 15236 MB successfully, after larger attempt(s) failed.
Error: Line 1 of updated_ref.csv has fewer tokens than expected.

My updated_ref.csv file looks like this
Name,chrpos
rs3793692,chr10_100008436
rs1983865,chr10_100016339
newrs150271620,chr10_100017402
rs10736126,chr10_100025924
chr10_100044928_A_G,chr10_100044928
rs7905091,chr10_100046078
rs6584182,chr10_100055644
chr10_100058710_C_T,chr10_100058710
rs4360594,chr10_10007547
rs927739,chr10_100077732
rs4917817,chr10_100085888

while my .bim looks like this
1 rs3094315 0 752566 G A
1 rs3131972 0 752721 A G
1 chr1770377 0 770377 T A
1 rs11240777 0 798959 A G
1 chr1834753 0 834753 C G
1 rs72890788 0 836924 A G
1 rs4970383 0 838555 A C
1 chr1840753 0 840753 G A

1 rs4475691 0 846808 A G

The reason is (likely) that plink is expecting a space separated file for the --update-name argument:
--update-name <filename> [new ID col. number] [old ID col.] [skip]

You've given it a comma-separated file, so it thinks there is only one column, hence the error. You also don't want column names.
Change the sep by:
sed 's/,/ /g' updated_ref.csv > updated_ref.txt

Or, edit the file in place, and keep a backup of the original named updated_ref.csv.bak :
sed -i.bak 's/,/ /g' updated_ref.csv > updated_ref.txt

Do you mean cat updated_ref.csv | sed -i 's/,/ /g' > updated_ref.txt
according to my data and it will convert my .csv into .txt?

See update to my post


text:
Precisely clipping bam file to bed coordinates
I have a bed file of coordinates and would like to clip anything outside these regions in my aligned bam files. However, trying programs such as bedtools intersect and samtools ampliconclip (on complementary bed file), they always let some reads slip through outside the ranges. I tried using hard-clip and both-ends options for ampliconclip aswell.
Example (top is clipped bam, middle is unclipped bam and bottom is bed file with regions I would like to keep):

Although most reads outside the designed amplicons are removed, there are still reads outside the regions in the Designed_pool1 file.
The problem is that the bed file contains amplicons and outside the amplicons are primers where I do not want to call variants. I cannot get the exact primer sequences because I don't know their length and the manufacturer does not state them, so I cannot use programs that require primer sequences.
Does anyone have some experience with this and know why the programs behave this way?
EDIT:
The exact commands I have used are:
bedtools intersect -a $BamFile -b $BedFile > $BamOut

and
samtools ampliconclip --hard-clip --both-ends -b $ComplBedFile $BamFile > $BamOut

as well as without --hard-clip and --both-ends options. All methods let some sequences stay out of the regions in bed files.
I am not completely familiar with the bedpe format, but I am not using paired end data. I cannot get primer start and stop positions, only amplicon positions.

Maybe you need to sort the bam file before trying to clip.  Did you try that already?

Yes, the bam files are sorted

I am curious whether your pileup visualization tool is respecting the clipping. It is possible that there is a subtlety of samtools ampliconclip that your visualization is not understanding, with regard to how it treats e.g. supplementary alignments or something.
Perhaps more importantly, I notice that the reads that are not being properly clipped in your visualization tend to have one end anchored in each of your BED regions.
Looking at the samtools ampliconclip docs, it seems to be assuming that reads will only be clipped from the ends. They will not be clipped if the region outside the BED is in the middle, if I read it correctly. I think that you may need to include an additional step in which you either trim reads (so that they are not long enough to bridge your regions) or you somehow postprocess your alignments to explicitly clip in the middle.
If true, this seems like a significant enough oversight of samtools that it may be worth logging an issue with the samtools maintenance team.
Also, not clear if it's relevant but reading the samtools ampliconclip docs:

adjustments to the left most mapping position (POS) will mean that coordinate sorted files will need resorting. In such cases the sorting order in the header is set to unknown. Clipping of reads results in template length (TLEN) being incorrect. This can be corrected by samtools fixmates. Any MD and NM aux tags will also be incorrect*, which can be fixed by samtools calmd.

For single-end reads, generically:
$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --fraction-ref 1 reads.bed roi.bed > answer.bed

The BED file answer.bed contains filtered reads which entirely overlap intervals in roi.bed. It contains fields required for backconversion to SAM/BAM.


text:
querying metadata (GDC Data Portal) in python3
I am trying to query metadata from the Genomic Data Commons (GDC) Data Portal and used the sample code from its website:
import requests
import json

cases_endpt = 'https://api.gdc.cancer.gov/cases'

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = ["submitter_id",
          "case_id",
          "primary_site",
          "disease_type",
          "diagnoses.vital_status"
         ]

fields = ','.join(fields)

params = {"fields": fields,
          "format": "TSV",
          "size": "100"
         }

response = requests.get(cases_endpt, params = params)

print(response.content)

However, I am not getting a tsv file in return. I am ultimately getting an FPKM instead. Any help is appreciated!

You can use str() to convert bytes to strings and print() with the proper separator.
import requests
from io import StringIO
import pandas as pd

cases_endpt = 'https://api.gdc.cancer.gov/cases'

fields = ["submitter_id",
          "case_id",
          "primary_site",
          "disease_type",
          "diagnoses.vital_status"
          ]

fields = ','.join(fields)

params = {"fields": fields,
          "format": "TSV",
          "size": "100"
          }

response = requests.get(cases_endpt, params=params)

content = str(response.content, 'utf-8')

print(content, sep='\t', end='\r')

If you want to convert to a data frame you can use StringIO() to load a data frame from an input string directly to a pandas data frame.
data = StringIO(content)
df = pd.read_csv(data, sep='\t')

The output:

case_id
disease_type
id
primary_site
submitter_id

0
375436b3-66ac-4d5e-b495-18a96d812a69
Adenomas and Adenocarcinomas
375436b3-66ac-4d5e-b495-18a96d812a69
Rectum
TCGA-F5-6810

1
38f93ca8-4971-4932-b592-0a17eaece2ad
Ductal and Lobular Neoplasms
38f93ca8-4971-4932-b592-0a17eaece2ad
Pancreas
C3L-02115

2
5ea1361a-050d-4fa0-be59-b530c456b6e3
Ductal and Lobular Neoplasms
5ea1361a-050d-4fa0-be59-b530c456b6e3
Pancreas
C3L-00928

3
690378f3-e38a-4b25-8489-05a3b4bd40b6
Adenomas and Adenocarcinomas
690378f3-e38a-4b25-8489-05a3b4bd40b6
Kidney
C3N-01175

4
6b4f086e-4ef9-4f30-874e-8a155b2e3bab
Squamous Cell Neoplasms
6b4f086e-4ef9-4f30-874e-8a155b2e3bab
Other and ill-defined sites
C3N-04276


text:
How to get KEGG pathway names from a list of 'ko' IDs?
I have a list of ko IDs from which I want to get the pathway names. How can I get that? I have tried this link, but it does not give me any result. What to do?
  ko00010   ko00030 ko00061 ko00120 ko00121 ko00190 ko00195 ko00230 ko00250 ko00253 ko00270 ko00300 ko00330 ko00400 ko00401 ko00410 ko00430 ko00450 ko00460 ko00471 ko00472 ko00500 ko00511 ko00520 ko00521 ko00523 ko00531 ko00550 ko00562 ko00604 ko00620 ko00625 ko00640

Thanks

You can use the KEGG API:
ko_list="ko00010 ko00030 ko00061 ko00120 ko00121
         ko00190 ko00195 ko00230 ko00250 ko00253
         ko00270 ko00300 ko00330 ko00400 ko00401
         ko00410 ko00430 ko00450 ko00460 ko00471
         ko00472 ko00500 ko00511 ko00520 ko00521
         ko00523 ko00531 ko00550 ko00562 ko00604
         ko00620 ko00625 ko00640"

for i in ${ko_list}
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${i/ko/map}")
    echo -e "${i}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done

Or if your list is too long you can change the input and use a file with a list of ids.
while read line
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${line/ko/map}")
    echo -e "${line}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done < ko_list.txt

The result will be something like:
 ko00010    Glycolysis / Gluconeogenesis
 ko00030    Pentose phosphate pathway
 ko00061    Fatty acid biosynthesis
 ko00120    Primary bile acid biosynthesis
 ko00121    Secondary bile acid biosynthesis

Thanks a lot Zorbax, for your answer. It works. Besides that, is there any web server or something that can be used? Just wanted to know.

Hi Zorbax, I have a similar query. Is there a similar way I can use to get KEGG Module names from a list of Module IDs (e.g. M00002)? Many thanks

I don't know anything about a web server. And yes, to your second question, I think you can use the API request for the module instead pathway list.

This doesn't work anymore

Yes, it does. Use https instead http


text:
How to deal with wildcards (N) in sequence logos
This is an open question that came to my mind recently, though it might be that it is not a common use case. The purpose is to find common patterns in a collection of sequences by representing the patterns as PSSMs. However, some of the bases have sequencing errors and are thus missing. Is there a way of integrating these missing bases in the PSSM calculation?
For a collection of sequences, deriving a sequence logo is relatively straightforward. For the case:
ACTAGCGT
ACGTACTC
ACGATCTC
ACGTCAGT

Taking the frequencies of each letter at each position would yield the following (pseudocount-corrected) Position Frequency Matrix:
      A       G     T       C
0  0.85  0.0375  0.05  0.0625
1  0.05  0.0375  0.05  0.8625
2  0.05  0.6375  0.25  0.0625
3  0.45  0.0375  0.45  0.0625
4  0.25  0.2375  0.25  0.2625
5  0.25  0.0375  0.05  0.6625
6  0.05  0.4375  0.45  0.0625
7  0.05  0.0375  0.45  0.4625

And by including background frequencies, we can derive the total Infomation Content for each position, and sum them up to compute the total Infomation Content of the matrix.
But what would happen if the analyzed sequences had experimented errors in the sequencing and some of the letters were missing?
ACNNGCGT
ACGTACTC
ACGATNTC
ACNTCAGT

Would you consider N a new letter of the alphabet and follow the same procedure? (counting ocurrences per position, calculating background frequencies, include it in the sequence logo...)
Would you ignore them and normalize each column's frequencies to the number of non-N bases?
Would you penalize  positions/sequences bearing N's when calculating the IC?
Other?

I wanted to start with the case of N since it's the most generic wildcard, and therefore could be easier to address; but if you have ideas on incorporating other more specific wildcards, feel free to add it :)

Since the Position Frequency Matrix is an estimate for the probability to find a particular letter at each position, I think the best way to deal with sequencing errors that result in N is to ignore that base for that position and normalize based on the number of non-N bases.
The Information Content at each position is just a function of the (estimated) probability distribution. So, once you have the 'best' estimate for the Position Frequency Matrix, I think you should use it directly with no penalty when calculating the IC.
If there is a position with a lot of Ns, though, you might get a biased estimate of the IC for that position - since using observed frequencies to estimate entropy/information is strongly biased when the number of observations is small.


text:
Understanding "Measurement of mRNA abundance using RNA-seq data" by Wagner et al 2012
I'm trying to understand why RPKM is not an appropriate way to normalize for RNA seq (I understand the general idea, but I'd like to gain a deeper understanding). So I'm reading the original paper by Wagner et al 2012 that suggests TPM as an alternative. However, I'm struggling to understand the following paragraph, and was wondering if anyone can shed some light on what the authors are trying to say:

The reason for the inconsistency of RPKM across samples arises from the normalization by the total number of reads. While rmc as well as qPCR results are ratios of transcript concentrations, the RPKM normalizes a proxy for transcript number by $r_g \times 10^3 / fl_g$ the number of sequencing reads in millions, $R / 10^6$. The latter, however, is not a measure of total transcript number. The relationship between $R$ and the total number of transcripts sampled depends on the size distribution of RNA transcripts, which can differ between samples. In a sample with, on average, longer transcripts the same number of reads represents fewer transcripts.

One part that I'm particularly confused by is: "The RPKM normalizes a proxy for transcript number by $r_g \times 10^3/fl_g$ the number of sequencing reads in millions, $R/10^6$." What is the "proxy" here?

That is why: https://www.biostars.org/p/9465851/#9465854.
All these naive per-million scaling methods may fail to correct for library composition differences.

Thank you for the reply. So do you mind elaborating on this idea of "library composition differences"? I keep encountering this idea when I try to understand  normalization but I just can't seem to understand it

Nicely illustrated here: https://www.youtube.com/watch?v=UFB993xufUU

Will check it out thank you. Another question: Do you know what the authors mean by "The RPKM normalizes a proxy for transcript number by $r_g \times 10^3/fl_g$ the number of sequencing reads in millions, $R/10^6$."
I understand that TPM is also not recommended to use, but I believe before I understand more complicated techniques I should build up my understanding using the simpler methods


text:
How to tell if our ligand-protein docking is good from AutoDock Vina's result
I have perform a ligand-protein docking using Autodock Vina.
The result of the docking looks like this:
WARNING: The search space volume > 27000 Angstrom^3 (See FAQ)
Detected 8 CPUs
Reading input ... done.
Setting up the scoring function ... done.
Analyzing the binding site ... done.
Using random seed: -1553787135
Performing search ... done.
Refining results ... done.

mode |   affinity | dist from best mode
     | (kcal/mol) | rmsd l.b.| rmsd u.b.
-----+------------+----------+----------
   1         -5.9      0.000      0.000.  Pose 1
   2         -5.7     22.945     25.492.  Pose 2
   3         -5.5      1.426      2.046.  Pose 3
   4         -5.5     23.669     25.616
   5         -5.4     25.783     29.152.  .....
   6         -5.3     21.146     23.357
   7         -5.2     20.323     22.545
   8         -5.2     23.864     26.064
   9         -5.1     23.422     26.585.  Pose 9

As far as I understand from these statistics Mode 1(Pose 1) is the best.
However when I actually visualize them in Pymol, Pose 1 has no hydrogen bonding at all
but Pose 2 has.
My question is how can we judge if which of those two Pose is the best to use?
Note in figure below Pose 2 has dashed line (Hydrogen bond).

It is best to contextualise the numbers. -1 kcal/mol is about the potential energy gained from a hydrogen bond —technically described in the r^6 part of the Lenard–Jones term, it is also the average collision energy of a water molecule at 37°C as that is RT ($\frac{k_b\cdot T}{N_A}$, wiki)under a Maxwell–Boltzmann distribution. A salt bridge –2 kcal/mol (Columbic force term).
So your scores are not very low, hence why you are counting two hydrogen bonds. Although you can also see a lovely sulfur-pi interaction which is –1 to –2 kcal/mol, so those bonds alone are probably making a –4 kcal/mol contribution, so I am guessing that some terms may be horrendous, such as repulsion forces etc.
A nice metric is doing a conversion to ligand efficiency, which weeds out affinities driven by size of the molecule... in this case most of the molecule is doing nothing.
Also, most programs have an accuracy of 1 kcal/mol or higher.
So one cannot say what ∆∆G is the best with the data at hand due to noise, but one can say that the ∆∆G is not low enough... Sorry.

Thanks. By  ∆∆G Do you mean  affinity (kcal/mol)?
Then what do you consider ideal score range for  ∆∆G. Since you said that my score above is not low enough?

I did and it generally is fine, but variations of definitions and calculation come into play. Gibbs free energy is ∆G (so one delta to start with), although some tools refer to ∆G_bind as ∆∆G, and generally it's intended as ∆∆G = ∆G_holo - (∆G_apo + ∆G_aq), where ∆G_aq is the free energy in aqueous solvent, but this is not always true as there may be fancy stuff correcting for ∆G_aq (which makes little sense for a high partition coefficient compound), Vina is one of them.

Higher than 10 kcal/mol... If you calculate the ligand efficiency it will be something between –0.1 and –0.2 kcal/mol/heavy. If your ligand is meant to dock with a high affinity, then I would revisit the template —a holo enzyme structure with the substrate removed is best, check the protonation of that aspartate etc.

Ah. I just noticed your box is giving a warning of being over 30 Å x 30 Å x 30 Å. Which is huge.

when you said 10 kcal/mol, how did you come up with that number?
Any reference/publication for that?

It's a ballpark figure. Playing devil's advocate, if you take the example of aspirin (an old small drug that's totally non-specific, that is taken at a 1,000 mg dose, that is charged, that is quickly broken by a liver cytochrome etc), it's binding to it's target is only -7 kcal/mol, which for 13 heavy atoms makes it's LE -0.5 kcal/mol, meaning that threshold doesn't always hold true (with massive caveats)


text:
Why isn’t ace2 included in SARS‑CoV‑2 PPI map?
I'm doing my master's thesis on drug repositioning for SARS‑CoV‑2 using AI more specifically network medicine,
I found only one human protein-protein interactions (PPIs) map from Gordon et al.
The manuscript doesn't include the ace2 protein which the SARS‑CoV‑2 spike protein directly binds with. AVE2 is a host cell surface receptor facilitating virus entry and thereafter replication.
Thus, I'm looking for PPI map that includes ACE2 protein but i can't find any
Secondly, is there a reason why it's not been included ?
PS: I have very basic knowledge of this domain

In fairness to the authors they do included it, viz. "captopril  modulates  ACE2  indirectly  via  its  direct  interaction  with  Angiotensin Converting  Enzyme,  ACE". So they hold results on ACE and therein have a claim for ACE2


text:
Sending chimera MD to remote computer
I am using Terminal to try to run Chimera MD simulations on a remote server/computer. I have been able to successfully connect using the command ssh username@detroit05 (and then entering the password). However, I want to be able to send a command to a remote computer to improve the computational time especially when running multiple MD simulations; right now, I am running MD simulations on the same computer. I was first thinking of a sequence of commands below:
Ctrl+Z
bg
ssh -y username@detroit05

I know that the Ctrl+Z command temporarily pauses a command, but I believe that the jobs command is first needed to list the current running jobs (that Ctrl+Z can be performed on). However, the running jobs are not listed when I type “jobs”. After using the “top” command to list all of the processes, I have verified that the Kill command works by typing these in succession:
Kill -STOP PID#
Kill -CONT PID#

Can you walk me through how to first have a listing of jobs running in Chimera and then sending it to a remote computer?


text:
Seurat - subsetting by genes expressed
Originally- I was looking if in at least one of these genes were expressed:
test1 <- subset(myData,
                  (KLRC1 > 0) | (KIR2DL1 > 0) |
                (KIR2DL3 > 0) | (KIR2DL4 > 0) |
                (KIR3DL1 > 0) | (KIR3DL2 > 0) | (KIR3DL3 > 0), slot = "data")

Now- I would like to subset if expressed in at least two of these genes- is this possible with the subset function?

With a little bit of workaround:
i) Add a new column to the data slot (only because your original subset() call  does so but it can be raw counts or any other data matrix in your Seurat object). rowSums() determines how many non-zero counts you have.
df$non_zero <- rowSums(cbind(df$gene1 > 0, df$gene2 > 0, df$gene3 > 0, ...)) > 1

Change df with the your data slot, should be sth like YourSeuratObject@...$data.
ii) Subset the newly modified data slot, basically treating the new column as another gene. Be careful not to use this column in preprocessing steps like normalization.
YourSeuratObject <- subset(YourSeuratObject, non_zero > 1), slot = "data")

Hi, Thanks so much for your quick response. I've tried many many combinations with my seurat object data slot -myData@assays$RNA@data$KIR2DL1 - but I dont know how to access the information via the gene name. Thanks again!!

Hi, Does anyone know how to access the gene name via the Seurat data sctructure as descibed in the answer above.

https://github.com/satijalab/seurat/wiki/Seurat#slots goes over how to interact with Seurat objects. Think of a Seurat object as a container and within it there are several assays (RNA, DNA, ...) and within each of these assays there are different (data) slots, for example one for raw counts, one for normalized counts... You might need sth like: GetAssayData(object = seurat_object, slot = 'data')[1:5, 1:5]. Check the manual for GetAssayData(), you can get a specific data slot from a specific assay.


text:
10X single cell sequencing issue
Can you please help me in a technical issue? We sent some samples to single cell sequencing The company says

Many of the reads were not assigned to cell-associated barcodes. This
could be caused by high levels of ambient RNA or by a significant
population of cells with a low RNA content, which the algorithm did
not call as cells. The latter case can be addressed by inspecting the
data to determine the appropriate cell count and using --force-cells.
This would happen if the cells aren’t transcriptionally active – but
as these are PBMCs this is unusual .

Have you ever heard about such a problem? Any solution in this case? sending more sample of each patient for sequencing or comparing single cell with matched bulk RNA-seq which one help?
10X say they would do deep sequencing (covering the cost) to get something for us but I don't know why my boss says we should compare single cell with matched bulk RNA-seq.

This is not a bioinformatics problem, it's a sample/library prep problem.

Are these single cells or single nuclei? Fresh or preserved in some way?

Agreed, but further bioinformatics might help better diagnose the issue, e.g. fastQC

@Chris_Rands Do we expect a different Fraction of reads in cell?

From what was described, mostly likely they ran cellranger and usually you can get a summary report like this
One of the metric would be Fraction Reads in Cells and I think this is what they are referring to this.
My suggestion would be to ask for the report from cellranger and check the metrics and the different plots. If what you have are indeed pbmcs, the reason is most likely not a lot of cells made it into droplets.
Submitting another sample will not help unless you can ensure the cells are of good enough quality, otherwise you will not get much cells back
Sequencing it deeper does not help as well because you are basically sequencing RNA of so called ambient DNA, and it does not give you useful information.

Answer from @devon-ryan, converted from comment:
Comparing with matched bulk RNA-seq isn't reasonably related to this, unless you want to ensure the pseudo-bulk is similarish to the actual bulk. Of course they'll be different, but I suppose if you had enough samples and enough background knowledge about what to expect maybe you could get an idea whether the pseudo-bulk way reasonable, which would hit at whether the scRNA-seq was OK. But realistically you'd need a lot of experience with these exact cells in bulk and scRNA-seq to do that. Otherwise, it's a waste of time.


text:
Querying metadata (GDC) using a filter
I have been trying to access the GDC metadata with python through the filter below. However, this code I have been using always yields a blank tsv file. Any help on fixing the code is appreciated!
import requests
import json

cases_endpt = 'https://api.gdc.cancer.gov/cases'

# The 'fields' parameter is passed as a comma-separated string of single names.
fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions"
    ]

    fields = ','.join(fields)

filters = {
    "op": "in",
    "content": {
        "field": "project_id",
        "value": ["TCGA-BRCA"]
    }
 }

params = {
    "filters": json.dumps(filters),
    "fields": fields,
    "format": "TSV",
    "size": "10000"
    }

response = requests.get(cases_endpt, params=params)

file = open("query10.tsv", "w")
file.write(response.text)
file.close()

I think you are using an outdated version of the API, you need to change the "field" with cases.project.project_id
import requests
import json

cases_endpt = "https://api.gdc.cancer.gov/cases"

fields = [
    "sample_id",
    "sample_uuid",
    "sample_type",
    "sample_type_id",
    "tissue_type",
    "tumor_code",
    "tumor_code_id",
    "oct_embedded",
    "shortest_dimension",
    "intermediate dimension",
    "longest dimension",
    "is_ffpe",
    "pathology_report_uuid",
    "tumor_descriptor",
    "current weight",
    "initial weight",
    "composition",
    "time_between_clamping_and_freezing",
    "time_between_excision_and_freezing",
    "days_to_sample_procurement",
    "freezing_method",
    "preservation_method",
    "days_to_collection",
    "portions",
]

fields = ",".join(fields)

filters = {
    "op": "in",
    "content": {"field": "cases.project.project_id",
                "value": ["TCGA-BRCA"]},
}

params = {"filters": json.dumps(filters),
          "format": "TSV",
          "size": "10000"}

response = requests.get(cases_endpt, params=params)
content = str(response.content, "utf-8")

with open("query10.tsv", "w") as file:
    file.write(content)

And if you use with open, you don't need to use close.


text:
How to create Beta Diversity distance boxplot
I'm trying to generate boxplot from Beta diversity boxplot between two conditions from a phyloseq object.
I have a Morisita-horn distance generated by vegdist fonction :
vegdist<-vegdist(t(otu_table(physeq_genus_3a)), method="horn")#matrice distance Morisita-horn
vegdist<-as.matrix(vegdist)

Now I have a matrix and I would like to somewhat create a mean distance boxplot of the Morisita-horn index winthin two condition like the graph below.
I wanted to do a simple boxplot with ggplot but beacause it is a matrix and not a dataframe, I can't do it.
Any help or suggestions would be more than welcome.

Could you please show a bit more of your data, or another test dataset? It would be greate if you could provide enough that we could use it to create a bar plot with error bars.

You can use as.data.frame() or as_tibble() to convert/coerce your matrix to a data frame or a tibble respectively (Technically a tibble is also a data frame). From there you should be fine with ggplot().

I would not use a "dynamite plot" (this is a box plot with error bars like the one in your question), that type of plots conceal the data! This post explains why and remedies.


text:
How to decide whether to assume 2 Gene Losses or 1 to find minimum cost trees for Gene Duplication/Expression events?
This paper introduces an approach to estimate a convergence from discord between phylogenetic species trees and gene trees, using Gene Duplication and Expression events.
A Gene Duplication (GD) event refers to when a bifurcation of the phylogenetic gene tree occurs creating two duplicate gene loci. The first image marks them in solid white diamonds.
A Gene Expression (GE) event refers to the loss or Regulatory Mutation (RM) of genes. The types are:

Gene Loss (marked by the cross in the first image) refers to when a locus is lost/deleted along the lineage of an extant species.
Gene Inactivation (marked by the unfilled downward caret in the first image) refers to a Regulatory Mutation where a gene is rendered inactive.
Gene Reactivation (marked by the unfilled upward caret in the first image) refers to a Regulatory Mutation where a previously inactivated gene is reactivated and sampled in extant species.

Now, in the event of a discord between a gene tree and a species tree, we can attempt to establish concordance through assumptions about hypothetical duplication and expression events. The gene tree with a minimum count of (GD + GE) is the most likely assumption for concordance.
In the first image, the right locus resulting from the second bifurcation of the gene tree is sampled in extant species C and D and undergoes a loss event for species E. I would assume that it would also undergo a loss event for species A and B but the paper merges those two loss events into one.
On the other hand, based on the second image from this paper, I've simulated the Gene Duplication/Expression events on the gene tree and found that if such loss events are merged into one, the Duplication/Expression (D/E in the table) count I calculate is lower than what is mentioned in the table.
I'm confused as to whether there is an appropriate approach to this. The author of the first paper mentions that such assumptions regarding Gene Duplication and Expression events are not at all representative of real Gene Duplication and Expression events (established through observed evidence) and serve as a mechanism to approach likely solutions to discordance between species trees and gene trees.
My question is, how do I arrive at the minimum cost tree?


text:
Unifying these files together
I have two .csv files
One for raw read counts of each gene and the other 3114 single cells for four patients which should be the row names of the first file
As single cells per patients seems duplicated I am getting error and don't know how to read them together
Any way this is the files
https://www.dropbox.com/s/slkhspz7lgdgph6/ascites.csv?dl=0
https://www.dropbox.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv?dl=0
I tried
a=read.csv("ascites.csv")
Warning message:
In doTryCatch(return(expr), name, parentenv, handler) :
  restarting interrupted promise evaluation
b=read.csv("ascitesRownames.csv")
View(b)
View(a)
rownames(a)=b

Error in .rowNamesDF<-(x, value = value) :    duplicate 'row.names'
are not allowed In addition: Warning message: non-unique value when
setting 'row.names': ‘logical(0)’

Any help in having a matrix in row names are 3114 single cells and genes in columns?

Can you please dput() the head of files so we don't have to download them from dropbox? Thanks.

There are a couple of problems:

As Ram RS said you get a data frame when you load that second CSV.  Since the information is stored in a row rather than a column you need a little more work to extract it too.
R will silently mangle your column names from the first data frame and will also interpret the second file as having column names (and then mangle those too) unless you give extra arguments to read.csv to prevent it.
The fact that you only have four identifiers (not 3114) in the second file means you can't use those as row names.  Your rows in a are cells, not patients, right?  You could make up something unique based on both patient and cell info but I'd suggest avoiding that.

First you could use something like this (going straight from Dropbox here) to read the files without losing information:
a <- read.csv("http://dl.dropboxusercontent.com/s/slkhspz7lgdgph6/ascites.csv", check.names = FALSE)
b <- read.csv("http://dl.dropboxusercontent.com/s/iq8lt7qrh7hdtq8/ascitesRownames.csv", header = FALSE)

Then, since you want to associate more information than can easily fit in a "wide-format" (matrix-like) data frame, and row and column names are so easily mangled by R, I'd suggest using a "long-format" data frame instead.  The package reshape2 is very helpful for this.  Here you can put the patient IDs as an additional column in your first data frame, explicitly label the cells, and then make a new data frame where each row is a cell and gene combination for a particular patient:
a <- cbind(patient = as.character(b[1, ]), cell_id = 1:nrow(a), a)
gene_counts <- reshape2::melt(
  a, id.vars = 1:2, variable.name = "gene", value.name = "count")

So you get something like this:
> head(gene_counts)
  patient cell_id gene count
1   7873M       1 A1BG     1
2   7873M       2 A1BG     0
3   7873M       3 A1BG     2
4   7873M       4 A1BG     0
5   7873M       5 A1BG     1
6   7873M       6 A1BG     0

(Note the gene column is a factor there.)
I'd be very careful with row and column names since R is so prone to altering them without warning, and in general with assuming things will match up in the right order like the rows of the one table and the columns of the other.  Hadley Wickham has a good summary of the problems with row names in particular:
https://adv-r.hadley.nz/vectors-chap.html#rownames

+1 agree with the comment about row names. I think data.table has the right idea in doing away with them completely and using keys instead.

read.csv produces a data.frame. You'll need rownames(a) <- b$col_name, not rownames(a) <- b.