text:
Finding common and unique data set by comparing two files based on their column and to split the columns multiple strings to print in output
I have  very large sizes tab-delimited .vcf files and want to match these two / or 3 files based on their position and print to a new .csv file
File structures:
File_1: tab-delimited file (.vcf) and its as column names as follows

(line number 3439) #CHROM     POS     ID       REF      ALT      QUAL    FILTER
 INFO
File 2: same as file_1 column names
(line number 3407) #CHROM     POS     ID      REF     ALT     QUAL    FILTER
 INFO  FORMAT  SAMPLE_1  SAMPLE_2 .....
In file_2, column 7 (INFO) contain many substrings like
AC=46,2;AF=0.958,0.042;AN=48;DP=269;ExcessHet=3.0103;FS=0.000;InbreedingCoeff=0.5411;MLEAC=92,4;MLEAF=1.00,0.083;QD=25.36;SOR=2.488  from these strings , have to print  only information of AC=, AF=, AN=, DP=
Desired output files to generate:

position matched in two of files: common_position_matched.csv

if I have more than 3 files also, the output should be one file and important thing, if a position column (1-POS) only matched only in 2 files and 1st file line should be NA NA NA
file_1  CHROM  POS    REF  ALT     file_2 #CHROM  POS      REF  ALT  INFO

1       22    10511521  C     T         1    22  10511521   C    T   AC=46,2 AF=0.958,0.042 AN=48 DP=269
2       22   10510544   G     A         2    22  10510544   G    A   AC=49,2 AF=0.958,0.042 AN=89 DP=536
3       22   10515068  AGAT,T AGAT,AT   3    22  10515068   AAA  AAAGG,A,GAA AC=100  AF=0.958,0.042 AN=62 DP=123
4       22   10515118  A G,   TAA       4    22  10515118   AG,  TAA AC=32   AF=0.958,0.042 AN=45 DP=500
5       22   10515118  AAAG   A         5    22  10515118   AATG A   AC=50   AF=0.958,0.042 AN=49 DP=129

note: while doing matching, not removing the duplicates, because in the same position there may be an addition or sometimes it may be deletion.

unique position of each file, in tab-delimited

output: File1_unique.csv and File2_unique.csv etc.
so far was able to read the file and match them according to position and print the output, but was not able to write efficient code
import pandas as pd
df1 = pd.read_csv("File1_3.vcf",sep='\t',usecols = ['POS']) ## Reading file1
df2 = pd.read_csv("file2_3.vcf", sep="\t", usecols = ['POS']) ## Reading file2
df3 = pd.concat([df1,df2], sort=True) ## Combining both the dataframes
df4 = df3.drop_duplicates(keep=False) ## Dropping the duplicates (intersect)
df4.to_csv("c3-UniquePosition_of_bothData.csv", sep="\t", index=False, header=True) ## Writing the unique to both
df1_Uni_file1_c3Posi = pd.merge(df4, df1, on='POS', how='inner') ## Identifying the unique position of File1
df2_Uni_File2_c3Posi = pd.merge(df4, df2, on='POS', how='inner') ## Identifying the unique position of File2
df_File1_File2_common_c3Posi = pd.merge(df1, df2, on='POS', how='inner') # Identifying the common chr-position of File1 and File2```

Program 2: (giving original file without editting)
import pandas as pd
df1= pd.read_csv("File1_22.vcf.gz", sep="\t", skiprows=3438, usecols = [0,1,2,3,4])
df2 = pd.read_csv("File2_22.vcf.gz", sep="\t", skiprows=3406, usecols = [0,1,3,4,7])
#writing the output files
#df1.to_csv("File1_c22.csv")
#df2.to_csv("File2_c22.csv")
#mergeing
df3 = pd.merge(df1, df2, on='POS', how='inner', sort=True)
df3.to_csv("common_position.csv", sep=",", index=False, header=True)
#df3 = pd.concat([df1,df2], axis=1).to_csv('check1.csv') # this command join multiple output to single output

Could any one give efficient python pandas script, to do this
Thanks All

I suggest you reformat your post to make it clearer. Right now I feel reluctant to even read it as it is extensive and not well formatted. maybe this is just me, but probably not.

Beyond reformatting, it would be great if you can clarify the question. What are the tab-delimited files? They appear to be lists of mutations, but can you confirm that? And what exactly do you mean by position matching? You provide many examples but it would be great if you can clarify in text exactly what those examples are meant to show.

@ATpo and alizbar,  I had reformated the question as much I can, could you people give me some solution

The code looks fine, a few things appear missing, assigning headers/indexing. It would be useful to see the output you are getting, because it a complex flatfile database and it will take me a while to look between all the columns and code

One way to do this is to use a dictionary. You can read the first file line by line, and for each line, use "chr:position" as a key, and the corresponding line as the value. If you want to handle duplicates, you can make the value a list and store multiple lines with duplicated "chr:position" as elements in that list.
After constructing this dictionary, read the second file line by line. For each line, set "chr:position" as the so-called name of that line. Then using this "chr:position", you can easily find the corresponding line in the dictionary. Append the line from the second file to that line from the first file.
After finishing reading the second file, you are ready to write the new file to output. To do so, iterate over the keys in the dictionary, and write the value to each key to a new file. Remember that each value in that dictionary is now a new line.
I choose to use a dictionary because they are basically a hash table, meaning that using the key you can instantly read that value in the dictionary, instead of iterating every element in a list to find that element. This is very handy when you have a large dictionary.

Well, after writing that long paragraph I realized that you are basically trying to merge two VCF files. Why not checking out vcf-merge from VCFtools.

This can be done entirely in pandas, there's bits missing. Its quite a complicated database problem to read

@Michael, I agree with you totally. Sometimes in R I use left_join and many other joining functions. I think Pandas has the same functionality. But sometimes I need to join millions of lines together (VCF files can be 700000000x1000), reading the data all at once seems unrealistic. I do not know much about database but may be this can be handled by some special tools?


text:
Edit all the fasta headers using awk
kindly your comment highly appreciated
have more than 4k header sequences look like:
>LTR_retrotransposon100_Gypsy?
cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

>LTR_retrotransposon10011_ERVK
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon10011_ERVL-MaLR
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>LTR_retrotransposon100_Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

I would like to use awk or sed to remove everything between the underscore, including the underscore and replace it with /
plus add sequential RNA followed by #.
and the target output could be:
>RNA1#LTR/Gypsy?
ctcagcagcactatttacaatagccaagacatggaaacaacctaaatgtcttatcaatag atgactggataaaggagctgtggtatatctatacaatggaataccatcagccataaaaaa gaataaaatattgccatttgcagcaacatggatggacctggagattatcattctaaggga agtaagccagaaagagaaagaaaaataccatatgatatcacttatatgtggaggtaaaaa aaaaaaaaaagacacaaatcaatttatttgcaaaacatacaTGGActttcagacatagaa

>RNA2#LTR/ERVL-MaLR
taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA3#LTR/ERVK
    taatttgaaaagatacatgcaacccaatgttcatagcagcattatttacaattgccaaga tatggaagcaatctaagtgtccatcaacagatgaatggagagagaagatgtggtatatat acaataaaatactattctgtcataaaaatgaataaaattctgccatttgcaacaacatgg atagacctggagggtattatgctatgtgaaataagtcagacagagaaagacatatactat

>RNA4#LTR/Copia
    cgtgcaccccaatgttcatagcagcactatttacaatagccaagacatggaaacaaccta aatgtcctttgacaggtgactggataaagaagctgtggtatatttatacaatggaatgct actcagccataaaaaagaataagataatgccatttacagcaacatggatagacctggaga atgttattctaagtgaagtaagccagaaagaaaaagaaaaataccatatgctatcactta

Is this in fasta format, that is, does the header start with >, like so: >LTR_retrotransposon10011_ERVL-MaLR? It should not matter for this post, but it may be an important detail for later steps using other bioinformatics tools.

yes > included, thank you i will edit the question

EDIT after OP's update:
Try this Perl one-liner:
perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file

Here, the command line flags are:
-e: tells the perl interpreter to look for the code inline, rather than in a file with the script,
-p: loop over the input one line at a time, execute the code in the one-liner for each iteration of the loop, and print the result into STDOUT.
The code has regex replacement (see perlre for more details):
s{>([^_]+)_.*_}{>RNA${i}#${1}/}:
replace an underscore, followed by zero or more occurrences (* quantifier) of any character (.), followed by an underscore, with slash (/).
Keep > and the original header part before the underscore, captured by parens into $1 variable.
Add string RNA${i}, where $i is the counter of the fasta headers. The counter is incremented by 1 when the substitution occurs, that is, when the next header occurs.
The variables are used like so in the substitution: ${i}, ${1} to prevent wrong interpolation.
I am also using s{}{} as delimiters for readability (as in most cases) instead of the default s///.
EXAMPLE:
echo ">LTR_retrotransposon100_Gypsy?\nACGT\n>LTR_retrotransposon10011_ERVK\nTCGA" | perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}'

Prints:
>RNA1#LTR/Gypsy?
ACGT
>RNA2#LTR/ERVK
TCGA

Dear thank you for your response, kindly check the updated question

Your own answer doesn't seem to do what you ask as it doesn't add the sequential identifier to each sequence. Maybe try this?
awk -v rec=0 '{
    if($0 ~ "^>"){
        rec++
        sub(">", ">RNA" rec "#", $0)
        sub("_.*_", "/", $0);
    }
    print $0
}' yourfile.fa

sed -e 's:^\(.*\)_.*_\(.*\)$:\1/\2:' < input > output

source: Stack Exchange

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

sort: votes


text:
Running tophat dockerfile in background
I need to use tophat2 to align some sequences, and I wish to use the docker container. These are large sequences so it will take a long time, plus I'm working on a university server so the chances of being kicked off the server are very high. How can I use nohup and & to run tophat2 through docker, if I need some user interaction?

Do you have root access? Docker generally requires root (as it normally is systemctl service and hard/janky to use otherwise)... Does the server not have a queuing system? Or is this a classic furtive job in the login node of an oversubscribed cluster? Also if an admin may kill a nohup job, they may kill your docker job.

Assuming your user is in the group permitted to use docker:
docker pull quay.io/biocontainers/tophat:2.1.1--py27_3
nohup docker run -v /some/path:/data quay.io/biocontainers/tophat:2.1.1--py27_3 tophat2 ... &

The only thing to remember is to mount everything you'll need into the docker container (/data in the example above). If you need to debug things then you can add the -it flag to docker run and use bash as the command.
As an aside, you might find using screen or tmux preferable to nohup ... &, since you can reconnect to sessions with them.

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


text:
Variant calling for a subset of genes using whole gneome sequencing data
I have few 100 raw fastq files from whole-genome sequencing data and I would like to map these files to a set of genes only (and not whole genome) so as to find SNP's associated with them. Can anyone tell me what snp analysis pipeline is best for targeted variation and not whole genome?
I have this in mind:
Mapping the reads with bwa and then using varscan
Hope to hear for more options

Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against the full reference and then extract the regions you are interested in, e.g. with the region option of samtools view:
samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]
Check its manual for details.
While there is basically nothign wrong with VarScan it is quite old and not maintained anymore. I suggest you give bcftools a try (check its manual). Other maintained alternatives can be (among others) freebayes, strelka2, VarDict or the GATK from the Broad institute.

vSNP can be used to target specific regions of a genome. Pipeline is a 2 step process:

create VCF files (via Freebayes)
then in your case run all 100 VCF files to create SNP table

A filter file can be supplied to step 2.  By filtering all position ranges other than the genes of interest, step 2 will output a SNP table (.xlsx) with only your genes of interest.  Also, table will include annotation if GBK is provide.


text:
How to visualize called narrowPeak files in UCSC Genome browser or IGV?
I have called peaks using MACS2. Then I got a narrowPeak file like this.
1   4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80
1   4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
1   4571554 4572026 testONE_peak_3  72  .   4.16005 7.21000 4.85200 282
1   4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
1   4784848 4785061 testONE_peak_5  31  .   2.24218 3.19160 1.21569 133

Then I tried to visualize it in UCSC Genome browser and I got following error.
Error File 'testONE_peaks.narrowPeak' - Unrecognized format line 1 of file: 1 4491713 4491913 testONE_peak_1  37  .   2.86003 3.72580 1.66717 80 (note: chrom names are case sensitive, e.g.: correct: 'chr1', incorrect: 'Chr1', incorrect: '1')

So I added chr to the chromosome column instead of 1 (Chromosome number).
chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

When I tried to load the file again I got following error.
Error File 'testONE_peaks_NEW2.narrowPeak' - Error line 1 of custom track: thickStart out of range (chromStart to chromEnd, or 0 if no CDS)

Then I added following line to the top and I could load the file to UCSC Genome browser.
track type=narrowPeak name="Somite narrowPeak"
chr1    4491713 4491913 testONE_peak_1  37  .   2.86003 3.7258  1.66717 80
chr1    4492782 4493153 testONE_peak_2  27  .   2.22933 2.71107 0.82442 242
chr1    4571554 4572026 testONE_peak_3  72  .   4.16005 7.21    4.852   282
chr1    4739486 4739879 testONE_peak_4  64  .   3.63128 6.42948 4.14394 275
chr1    4784848 4785061 testONE_peak_5  31  .   2.24218 3.1916  1.21569 133

But the problem is I can not see anything in the Genome browser. Only the name is there. But nothing else. Do you know what I should do?

Please show a screenshot.

I am struggling with the same problem, where did you add this line of code?
awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > fixed_your_peaks.narrowPeaks

You can run this command in a command line or terminal

I think you might have changed the separator (or at least have some kind of inconsistency from the required format) for your file. Note that peak output files from MACS2 are variants of BED files.
It seems you need to have a tab separator for this type of file.
I copied your example file and then ran awk -v OFS="\t" '$1=$1' your_peaks.narrowPeaks > fixed_your_peaks.narrowPeaks on it. This produced a result that I was able to see in both the UCSC genome browser and IGV.

Strangely if you simply add ".bed" to your file it should work (i.e. testONE_peaks_NEW2.narrowPeak.bed)

Can we do this in R?

Oh. Thank you very much.


text:
Index a BAM file using pysam
(How) can you index a BAM file using pysam?
When I tried the intuitive pysam.index I got:
import pysam
my_bam = pysam.AlignmentFile("regular_bwamem_mapping.bam", "rb")
pysam.index(my_bam)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-5-a0b3a04ecb2c> in <module>
----> 1 pysam.index(L_bam)

~/miniconda3/lib/python3.7/site-packages/pysam/utils.py in __call__(self, *args, **kwargs)
     59             args,
     60             catch_stdout=kwargs.get("catch_stdout", True),
---> 61             save_stdout=kwargs.get("save_stdout", None))
     62
     63         if kwargs.get("split_lines", False):

pysam/libcutils.pyx in pysam.libcutils._pysam_dispatch()

TypeError: object of type 'pysam.libcalignmentfile.AlignmentFile' has no len()

disclaimer: when I post a question and answer, it's because I figure the answer while I write the question.

I have tried same command...
error
[E::hts_idx_push] Unsorted positions on sequence #8: 3411049 followed by 3378665
[E::sam_index] Read 'A00404:40:HFGYNDSXX:4:1258:24623:27633' with ref_name='chr8', ref_length=145138636, flags=1107, pos=3378665 cannot be indexed
Traceback (most recent call last): File "Insertx.py", line 283, in  main() File "Insertx.py", line 263, in main pysam.index(str(''.join( data))) File "/usr/local/lib/python3.8/dist-packages/pysam/utils.py", line 69, in call raise SamtoolsError(
pysam.utils.SamtoolsError: 'samtools returned with error 1: stdout=, stderr=samt

@anju there must have been something that gone wrong with sorting the file. Otherwise it would not complain about it.

Oh you silly sausage, pysam.index takes a bam file name, not a python object.
import pysam
pysam.index("regular_bwamem_mapping.bam")

will index your .bam file.

I should print "Oh you silly sausage" and tape it above my monitor.


text:
Why do these NCBI representative genomes for ape species have no Y chromosome?
I recently downloaded a genome from NCBI for chimpanzee and was surprised to see no Y chromosome.  I added in one from a separate accession and carried on, but then today went to do the same for gorilla and saw the same thing, no Y.

For Gorilla gorilla (western gorilla): https://www.ncbi.nlm.nih.gov/genome/2156
For Pan troglodytes (chimpanzee): https://www.ncbi.nlm.nih.gov/genome/202

Why would this be?  I can see a separate entry for a Y chromosome for that same exact chimp used for the representative genome (Clint) here, but I don't see why it's not in the representative genome listing itself:
https://www.ncbi.nlm.nih.gov/nuccore/DP000054.3

tl;dr: technical difficulties, or sex
For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome.
For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that does have a Y chromosome sequence assembly that is directly linked, they just haven't updated the main reference sequence yet.

To directly address the question, the Y chromosome is very hard to assemble. The easy part of it to assemble is already on the X chromosome anyways (the pseudoautosomal region), and the hard part of it is often a mess of repeats, especially in mammals. Many groups will decide to just assemble the X as well as they can and call it good, due to the difficulties associated with assembling the Y. If they are working with a male individual the Y sequence is still there, it just may not be anchored to a chromosome sequence, so it's only in the unscaffolded/unanchored assembly ("master WGS" I believe, usually). It is also probably highly fragmented (it isn't a whole chromosome of sequence stitched together, but rather a lot of little pieces of sequence).
In some cases they will focus all sequencing efforts on female individuals specifically to make their assembly job easier, because in females the X is at the same copy number as the autosomes and doesn't have interference from homologous regions on the Y.


text:
Detecting whether genome is genetically engineered?
Is there a methodology for detecting whether a genome has been tampered with by a biologist?  For instance, say there has been an ebola outbreak such as in "The Andromeda Strain".  Is there a way a researcher can tell by just looking at the ebola genetic data and determine it is a genetically engineered strain?
Some ideas I have:

organism is abnormally close to unrelated organisms, i.e. some of the genetic code from the other organisms has been inserted into the first
subsequences are very different from the rest of the genome.  If there is a consistent rate of mutation applied to the entire genome, then a subsequence that is abnormally noise free could be inserted.
combination of 1 and 2: finding clear instances of unrelated organism genomes inserted into the genome being examined.

Additionally, are there any publicly accessible reference genomes that are known to be genetically engineered, that I could use as examples to test such ideas?
UPDATE: the wikipedia article on detecting genetically modified organism is not quite what I am looking for, since it claims current methods cannot detect unknown genetic modification.
https://en.wikipedia.org/wiki/Detection_of_genetically_modified_organisms

Currently, it is highly unlikely that the presence of unexpected or
  even unknown GMOs will be detected, since either the DNA sequence of
  the transgene or its product, the protein, must be known for
  detection.

Researcher have generated transgenic organisms for over a century. So yes there are many "modified genomes" available. You could also take any strain of interest, modify a few nucleotides, add some indels, or even an entire gene and analyze that. For step 3 for example, you could take humanized mouse or fly models used in disease biology...

@MichaelG. since you reference conspiracy theories, probably in regard to the coronavirus outbreak: for the record, I do not think the Wuhan coronavirus is genetically engineered.  But, the conspiracy theories made me wonder if there is a reliable way to detect whether genetically engineering has occurred in genetic data.

For the record, there is no guaranteed way to do this, it'd be completely impossible. We genetically engineer things all the time by just changing a few bases to things that already naturally occur in the species (but not at the same time in the wild) to speed up animal/crop breeding, as a simple example.

@DevonRyan can we at least guarantee true positives with high reliability?  It seems to be the case, e.g. if we see a massive segment of DNA from a totally different species spliced into the DNA sequence we are examining.  I'm asking if there is a good general methodology for at least identifying true positives reliably, even if we cannot identify true negatives reliably.

Sure, transgenes will make things obvious, but they're not generally present when we genetically modify something.

The definition of a GMO for example UK GMO regulatory bodies (HSE) is strictly the artificial introduction of foreign genetic material using targeted genetic engineering, such as splicing,  genome integration etc... Attenuation for example by serial passage the basis of the yellow fever vaccine, does not count as a GMO.
The simple answer is look at phylogenetic trees between the GMO strain against the total distribution of isolates. If large regions have been spliced the 'silent mutation rate' looks very different,.because the GMO parent is from a very different environment to the recipient.  The switch in replication rate affects the numbers of 'silent mutation', which is fast in RNA viruses.
If the GMO is a point mutation, ie limited amino acid change then I've actually done this with dengue virus and it is fairly easy because the novel mutation make the GMO 'stick out'.of the phylogeny against all natural mutation patterns. At this point the outlier is tracked back to Genbank where the relavant publication defines the engineered virus.
If you are inferring a conspiracy theory that is a different process and a different discussion, particularly if it is Ebolavirus which is a Category 4 pathogen. The investigation is again within the phylogeny, but in this case the answer is always the isolate is a natural zoonose from a sylvatic (wild) transmission cycle which has spilled over into humans.

Is there any way to do this if all you have is the particular strain and no database to compare it against?

No there is no internal signature. You need to compare the putative engineered (virus) against all other isolates, wild-type viruses.

Hmm, but we can detect plagiarism in this way.  For example, let's say I have a long string of text and Lorem ipsum dolor sit amet, consectetur adipiscing elit.  :)  It seems like we could determine the characteristics of much of the sequence and then detect deviations from these characteristics, at least if the deviations are large enough.

Plagiarism is detected by reference to prior publication, same idea with DNA/RNA

I'm referring to stylometry: https://en.wikipedia.org/wiki/Plagiarism_detection#Stylometry.  If something like this is not done in bioinformatics, seems like a good research avenue.

If you need to engineer a trait, then it has not evolved, so it probably does not confer a selectional advantage to the organism.
So a fingerprint of an engineered organism is that the genetic difference between it and its relatives does not have a selectional explanation.

Most of the things in nature don't have a selectional explanation.


text:
Why the gene symbol aliases have significantly different sequences?
"Aliases" or "Synonyms" should represent the same gene but with different names. But when I try to find the sequence of one alias and try to match with another alias they are significantly different.
During bioinformatics analysis, I need to convert NAT5 alias ID to its gene symbol ID (NAA20, NAA50). When I tried to match the sequences of NAA20 and NAA50, they are very different but still NAT5 is an alias of both NAA20 and NAA50.
I have to convert one id to another id (ENSEMBL). These gene symbol and alias is so confusing.
I am looking for NAT5 and when I search for NAT5 in ENSEMBL, it shows both NAA20 and NAA50.
From UniProt the canonical sequence of both the proteins:
https://www.uniprot.org/uniprot/P61599#sequences
https://www.uniprot.org/uniprot/Q9GZZ1#sequences
I am confused which one to chose?

At Ensembl, we categorise synonyms as anything that a gene might also be known as. This includes older names for them, since those names will be in the literature, including where a gene has been split in two.

Answer from @devon-ryan, converted from comment:
Both, that's why NAT5 is no longer a human gene symbol.
NAT5 does not equal NAA20, it's an out-dated name for it. Names change over time as people realize that there are more genes for something than originally thought.


text:
How to get the GO information for all the human genes?
I would like to retrieve a list of all the human genes and their proteins’ functional classification.
The Gene type attribute on ensembl.org is not specific enough, because it assign the label coding protein to all the proteins for which I am trying to classify. Also the KEGG Pathway and Enzyme ID option doesn't seem to help.
I have translated all the gene names to their GO representation, but now I am not sure how to query for their cellular function, biological process, and cellular component.
Is there such a platform that can retrieve this info from PANTHER?
Sadly, the reference proteome of humans doesn't contain a functional classification either:
 zcat UP000005640_9606.fasta.gz | grep -i mYOSIN-2
>sp|Q9UKX2|MYH2_HUMAN Myosin-2 OS=Homo sapiens OX=9606 GN=MYH2 PE=1 SV=1

# >db|UniqueIdentifier|EntryName ProteinName OS=OrganismName OX=OrganismIdentifier [GN=GeneName ]PE=ProteinExistence SV=SequenceVersion

Is there any way to use bgee, QuickGO, GONet, or Uniport to cluster genes using the customize columns?
EDIT:
I think that getting the most specific (GO term name) and least specific (GO domain) gene's GO values can be retrieved from ensmebl:

If you really need ALL genes then you won't find such a database, because we don't know all of their functions.

Have you looked at GO annotations? They are an hierarchical annotation for gene function. For enzymes the EC numbers are also an option. PFam is for structural domains, which generally align with function.

GO terms are used routinely to cluster genes for enrichment analysis. There are three categories, location, function and process, and a gene/protein will have one or more descendants (it's a hierarchical ontology). There are online tools that classify genes, such as GOrilla or Reactome. Programmatically you can use say Goatools in Python. There is a R package (don't recall name). And I think a Julia.

@MatteoFerla thanks. Let me know if you have a working example. I found GOATOOLS: https://www.nature.com/articles/s41598-018-28948-z


text:
How can I find genes of fundamental metabolic processes for an organism?
I want to find genes of fundamental metabolic processes in the Phytophthora cinnamomi genome. I'm using the genome sequence deposited in the NCBI and I'm searching ORF by ORF, using the ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/). Next, I use SMARTBLAST (https://blast.ncbi.nlm.nih.gov/smartblast/?LINK_LOC=BlastHomeLink) to find similar organisms.
How do I know that the gene that I'm searching expresses proteins of fundamental metabolic processes, and if the organism is similar to mine?
Can the best hits in SMARTBLAST tell me something about that?
Thanks

Can you please not SHOUT?

Should be a bit quiter now ;)


text:
Why my bim file doesn't match to my ped file as the Plink documentation suggests?
Plink documentation about .bim says the columns 5 and 6 are the Allele 1 and Allele 2 respectively.
The documentation about the .ped in turn, says "...The first six fields are the same as those in a .fam file. The seventh and eighth fields are allele calls for the first variant in the .map file ('0' = no call); the 9th and 10th are allele calls for the second variant; and so on.."
I've simulated a dataset as per:
$ plink1.9 --simulate ds1.sim acgt --out ds1
...
$ head ds1.bim
1       null_0  0       1       C       T
1       null_1  0       2       A       C
1       null_2  0       3       C       A
1       null_3  0       4       T       A
1       null_4  0       5       G       C
1       null_5  0       6       T       C
1       null_6  0       7       T       C
1       null_7  0       8       G       C
1       null_8  0       9       C       G
1       null_9  0       10      T       A

Than I converted it to plain PED and MAP:
$ plink1.9 --bfile ds1 --recode --out ds1
...
$ head ds1.ped > ds1.head.ped
$ vim ds1.head.ped
per0 per0 0 0 2 2 T T A C A A T A C C T C T C C C....

Given the data above, why my .bim file shows "C and T" for the first variant while my .ped file shows "T and T"?
Apparently they don't match as per described in the documentation (or more likely, I'm misinterpreting it)...
Also, how do those 2 files related to each other? E.g.: How do I know that my Sample 1 has the SNP's null_3, null_4 or null_5, etc?
I need to clarify those doubt because I want to fit a LightGBM model with the Genomic Data (can be plink or vcf). Not sure yet how to do. I was trying to convert the plink files to a CSV-style somehow, to be able to use in my model. Not sure if this is the way to go though...
Thanks!

Responded in the plink2-users group (https://groups.google.com/forum/#!topic/plink2-users/zgJxdXxvdLo ).


text:
Q: What analyses can I perform with a completely phased genome assembly?
Phasing a genome is the process of determining which variants lie on which copy of each chromosome. For example, position 5,430,500 might have an A on the paternal chromosome, but a T base on the maternal chromosome.
A haplotype block is a section of the genome for which we know which haplotype has what variants. A completely phased genome has phase blocks that are the full length of each chromosome. Every variant on every chromosome can be assigned specifically to the maternal to or the paternal haplotype.
My question is, what can we learn from having completely phased genomes that we cannot from an unphased genome, or from a phased genome with short haplotype blocks?

Some population genetic inference methods require the data to be phased because they look for long stretches of matching alleles between individuals, rather than just matching genotypes. Here, 'genotype' data is un-phased, and 'allele' data is phased.
This is because there is a degree of ambiguity when it comes to matching heterozygous positions - if the data is un-phased, you don't know if 2 individuals share mutations on the same chromosome (here, I mean different copies of the same chromosome - like either copy of chromosome 1, for example). This is important, because the fact that the mutation is shared on the same chromosome means that they are much more likely to share recent common ancestry with one another.
Let's say we have genotype data for two individuals (where the number corresponds to the number of non-reference alleles present):
Ind 1

0 ----- 1 ----- 2 ----- 0

Ind 2

0 ----- 1 ----- 2 ----- 0

If were looking at unphased genotype data only, we know the alleles match exactly in the 1st, 3rd and 4th positions, because they are homozgyous (either 0 or 2 copies of the non-referencee allele), but because the second position is heterozygous, we can't be sure which chromosome the non-reference allele is on. This is because the haplotype configurations for individual 1 could be either
0 ----- 0 ----- 1 ----- 0 Chromosome1a
0 ----- 1 ----- 1 ----- 0 Chromosome1b

or
0 ----- 1 ----- 1 ----- 0 Chromosome1a
0 ----- 0 ----- 1 ----- 0 Chromosome1b

i.e., we don't know which chromosome the non-reference allele at position 2 is on.
If we then phased the data, and found that the non-reference allele was on chromosome 1a for Individual 1, and chromosome 1b for Individual 2 there is less evidence of recent shared ancestry than otherwise. I think this is more or less the approach that 23andme takes when they are inferring your relationship to other individuals in the dataset. After writing this, I realised there is a much better description of what I just said on the 23andme website!
For example, the above applies to the ChromoPainter/fineSTRUCTURE algorithms. This means that they are able to detect more subtle genetic differences between individuals than if you were to use non-phased data.
This is a slightly niche example, and there occasions where you need phased data, but this is the one I am most familiar with :). Off the top of my head, methods such as Multiply Sequential Markovian Coalescent (MSMC), RELATE, and Extended Haplotype Homozygosity all require the data to be phased.

A situation where having phasing data is useful is with allele-specific protein binding.
Some individuals will have a heterozygous mutation within the binding site of a particular protein.
If you perform a standard chromatin immunoprecipitation (ChIP) sequencing assay, you'll find sequencing reads around the protein's binding site, but you won't know which allele they'll come from since sequencing is an ensemble measurement.
However, if you have phasing information, you can measure whether this specific mutation has an effect on that protein's binding, and get some insight into the function of that mutation in the context of that individual's cells.


text:
How to calculate Centimorgan lengths from run lengths in one-to-one dna profile comparison
For easier management of a database of dna profiles, I'm creating a tool to perform one-to-one matching on two profiles with snp genotype information.
Basically, something that works similar to GEDMatch's one-to-one matching feature.
Determining the length of runs of half or fully matched snp's, including some allowance for a genotyping error rate, is the easy part.
My question is about converting an snp run length to centimorgans.
I figure that the length needs to be adjusted by the likelihood of a change in the particular region of the chromosome.
Since my dna profiles and the positions of snps are based on Genome Reference Consortium Human Build 37, I figure that his information could be in the reference build, directly or derivable in some way.
Looking at the GRCh37 download files and the definitions of their formats, I'm unable to find such information.
Is this information present in GRCh37, or can it be downloaded from another reference source, and from which field in the dataset can it be retrieved or derived?

Try having a look at this link:
http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh37.map.zip
If you unzip the file, there is a file for each chromosome (1:22) and then 3 for the x chromosome. The first column of each file details which chromosome the file corresponds to, the 3rd column gives you the cumulative genetic distance in cM and the 4th column gives you the genomic coordinate (position) in bp.
You can then calculate the distance in cM between two positions from this data by subtracting the cM positions between the two positions you are interested in.


text:
How to initiate RNA-Seq analysis of TCGA files?
I want to compare RNA-Seq datasets obtained from the TCGA to investigate how my gene of interest is implicated in different types of cancer.
I'm trying to download the data from the GDC Data Portal (https://portal.gdc.cancer.gov), but it is not user-friendly.
Can somebody give me some initial directions? I found a bioconductor package that facilitates access to the TCGA (https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html), but first of all, I would like to understand how the files are organized, the nomenclature and how to inspect them.
thank you very much, Fabiano

When it comes to user-friendliness, I found the FireBrowse portal (http://firebrowse.org/) much easier to navigate than the GDC one.

thank you very much, vuzun. In my opinion, firebrowse is interesting to have some preliminary results, but for a deeper analysis and a high-quality publication TCGAbiolinks is still the best option.

The GDC API user guide (https://docs.gdc.cancer.gov/API/PDF/API_UG.pdf) provides some guidance on the structure of GDC data. So for getting RNA-Seq data from GDC, you would provide this filter to get the 551 text files with workflow_type set as "HTSeq - FPKM-UQ" and "TCGA-LUSC" for the project. There are other workflow types as well:
cases.project.program.name in ["TCGA"]
and cases.project.project_id in ["TCGA-LUSC"]
and files.access in ["open"]
and files.analysis.workflow_type in ["HTSeq - FPKM-UQ"]
and files.data_format in ["txt"]
and files.data_type in ["Gene Expression Quantification"]
and files.experimental_strategy in ["RNA-Seq"]

(You don’t need all those filters to get the 551 files; I included them to show what is available.) The txt files are returned as gzipped tar files that hold text files with rows of the form:
ENSG00000158486.12  5341.21601153
The above filters provide all the FPKM-UQ files for TCGA-LUSC. If you wanted to just get the data for normal samples as controls, that would require filtering on sample type as well, adding in the additional filter:
and cases.samples.sample_type in [“solid tissue normal”]

That will return 49 files. If you instead wished to get only the tumor files, the additional filter would instead be:
and cases.samples.sample_type in [“primary tumor”]

and that returns 502 files. Note that there are several different sample types for tumor and normal tissues; those just happen to be the ones used here.
(Note that if you are using the GDC Repository portal, adding that additional filter requires you to select “Add a Case/Biospecimen Filter” on the “Cases” filter tab.)
Since GDC is part of the NCI Cancer Research Data Commons (CRDC), instead of downloading the data, you can actually work with it in the cloud using the various CRDC Cloud Resources such as the ISB-CGC: https://www.isb-cgc.org/. In fact, all the RNA-Seq data for TCGA is in an ISB-CGC Google BigQuery table: (https://console.cloud.google.com/bigquery?p=isb-cgc&d=TCGA_hg38_data_v0&t=RNAseq_Gene_Expression&page=table)
And the above count data can be pulled using the query:
SELECT HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC"

If you want to discriminate between normal and tumor expression levels, the ISB-CGC BQ tables make it pretty easy to tag the expression values with the type of tissue involved. The sample type can be assigned by joining the previous query with the biospecimen table:
WITH bio as (
SELECT sample_barcode, sample_type_name
FROM `isb-cgc.TCGA_bioclin_v0.Biospecimen`
),
expr as (SELECT sample_barcode, HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC")
SELECT bio.sample_type_name, expr.HTSeq__FPKM_UQ
FROM expr
JOIN bio
ON bio.sample_barcode = expr.sample_barcode

The GDC has listed some of the available sample types on this page: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
Also, BigQuery has a user interface that makes exploring data sets and writing queries easy, but you need to log in with a Google ID, so the ISB-CGC has a BigQuery Table Search (https://isb-cgc.appspot.com/bq_meta_search/) that allows you to explore the available tables without logging in. If you want to run queries, Google offers a free tier that you can use to try things out and run small jobs (https://cloud.google.com/free/) and ISB-CGC offers free credits too. ISB-CGC also has examples of working in R with BigQuery within the Community Notebook Repository (https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowTos.html). A few useful notebooks would be 'How do I create cohorts of patients?', 'How do I plot a heatmap using data in BigQuery?', and 'How do I get started fast?'.

thank you very much for your great feedback, Lauren. It seems an excellent option to work in the cloud. The only problem is that I only know how to program in R.

and a little in Bash..

Yeah, writing SQL can take some getting used to, but you can do quite a bit of stuff with simple queries. I added some more information to the answer for working in R with BigQuery.

many thanks, Lauren. great answers!

Lauren, I' m having difficulties to find the healthy controls. As far as I know, the TCGA contains tumor samples and healthy tissue samples (tissue adjacent to the tumor that serves as a normal control). I need the normal controls to know if my gene is up regulated or down regulated. Where these controls are located?

The tissue type labels are in the biospecimen table. I updated the answer to include a query for joining the biospecimen table with the RNA sequencing table. Also, I changed the original project ID from TCGA-OV to TCGA-LUSC in the answer, since it turns out that TCGA-OV has no RNA-Seq data for normal samples.

I find UCSC Xena to be well-organized and easy to use: https://xenabrowser.net/datapages/
Plain text files that are clearly labeled for all TCGA data (that is not protected).


text:
Using cellranger for non-10x data
I'm trying to use cellranger to mkfastq and then count and aggregate single-cell data. This data was not generated using 10x Genomics RNASeq.
When I run count, cellranger fails to auto id the chemistry and crashes out. Is it possible to wrangle the data into the cellranger method or should I look elsewhere. I'd ideally like to use the same method for all my single-cell data but not all of it is 10x.

Hi Yuriy, could you tell us what how you generated the RNAseq data? That might be important. Also what species/problem you are working with?

How do the fastq files that belong to a given sample look like, I mean how are they named? Can you share the ls output of a given sample directory?

Perhaps STARsolo is worth looking into.

Hi, It is Smart-Seq2 data that I am hoping to overlap with some 10x data. Regarding the fastq I had to use the mkfastq command. For each sample, I have _R1 _R2 _I1 _I2 files. The naming consists of the sample name, the clone, the lane and the read.

Possible?  Sure.  Your biggest problem is that 10x is expecting barcodes from its whitelist.  I used to alter the python scripts to accept a 'new' kind of chemistry where the barcodes and umi lengths matched what I had, but I found that my modifications stopped working when the version updated,  and I couldn't figure it out, so I started modifying fastqs instead.

Do you mind sharing the altered code? And how do you modify the fastqs?

My modifications didn't work, and I strongly recommend not changing code like that unless you understand what you are doing.  As to changing the fastqs, you have to understand the format you've got, and what format you want.  I use umi_tools to identify the barcodes and their one-offs, and use that list to make replacements.


text:
Permanent error with maftools
I want to remove variant including in copy number changed area by maftools but I am permanently get this error
> tcga.ab.3009.het = inferHeterogeneity(maf = laml_res, tsb = tsb, segFile = seg, vafCol = 'T_VAF')
CN data for following samples not found. Ignoring them ..
 [1] "LP6008337-DNA_A07" "LP6008460-DNA_F02" "LP6008460-DNA_G03" "LP6008460-DNA_D01" "LP6008202-DNA_B03" "LP6008460-DNA_A04"
 [7] "LP6005935-DNA_G02" "LP6008334-DNA_B02" "LP6008269-DNA_F06" "LP6008336-DNA_F02" "LP6005334-DNA_H01" "LP6005500-DNA_D03"
[13] "LP6008334-DNA_A04" "LP6008336-DNA_G01" "LP6008336-DNA_H02" "LP6007600"
Processing LP6008334-DNA_D02..
Error in data.table::foverlaps(x = tsb.dat, y = seg, by.x = c("Chromosome",  :
  The last two columns in by.x should correspond to the 'start' and 'end' intervals in data.table 'x' and must be integer/numeric type.

> head(dat.tsb[1:3,1:6])
   Tumor_Sample_Barcode Chromosome Start_Position End_Position Reference_Allele Tumor_Seq_Allele2
1:    LP6008334-DNA_A03       chr5       70858261     70858261                -                 A
2:    LP6008334-DNA_A03      chr10      115987869    115987869                A                 -
3:    LP6008334-DNA_A03      chr16       29859285     29859287              ACC                 -
>

How can I do this?

What are the types of Start_Position and End_Position columns? You can get that info with str(your_table).

Thank you they are integer (start and end of genomic range where the mutation supposed to be happened)


text:
RepeatMasker annotation
I collected my own denovo library using many different tools Pipeline (structure-based ) and I updated all the headers corresponding to RepeatMasker format, to obtain Repeat coverage in the genome and .tbl annotation.
when I run repeatmasker the annotation results assume to be Homosapien, where Homo annotation table ( .tbl) is limited as shown below Table (1)
Actually my denovo repeat library consists of more than what is in the table above, example, rolling circles, Copia, Gypsy etc..
RepeatMasker -Lib myOwnLibrary.fa Mygenome.fa

Table (1)

My question is how can I run RepeatMasker with my own library using -lib Flag and result like the table (2) below,
the table below assume to be vertebrates genome, this is my aim.
RepeatMasker -Species Vertebrates Mygenome.fa

Table (2)

Regards

Hello, a legit question. Could you also add the command you used to run the repeatMasker? Also, what is the second table? is it another repeatMasker run you did?

yes. it is another run, I updated the question

After long years of just conceptual knowledge of what the software does, I was dealing with RepearMasker this week. So, I am not entirely sure I am correct, but this is how I understand it.
The two flags you explored serve for different scenarios:

-species - if the taxa/species is present in the repeat database (repBase or Dfam)
-lib - if you want to give repeatModeler something different.

In your case it seems like your library is completely missing all LINEs, SINEs and a few others. What you can do is to get from the database all repetitions of the corresponding taxa and catinate it with the repeats you are interested in (the one you supplied in the -lib argument in the first command).
util/queryRepeatDatabase.pl -species Vertebrates > flavoured_vertebrate_repetitions.lib
cat myOwnLibrary.fa >> flavoured_vertebrate_repetitions.lib
RepeatMasker -lib flavoured_vertebrate_repetitions.lib Mygenome.fa

This should produce annotation similar to the bottom case, but also using the repetitions you manually added.
P.S. util/queryRepeatDatabase.pl is a script to query the database of repeats and it is shipped together with repeatModeler.
-- edit --
The TE classification in the summary file of repeatMasker is based on headers of the repetitive elements in the fasta file you specify in to -lib argument. So if you want to get the second-table-like, you need to somehow reclassify the repeat library to contan repBase-like classification of TEs. The explanation of lib argument from the manual:

The recommended format for IDs in a custom library is:
>repeatname#class/subclass or simply
>repeatname#class
In this format, the data will be processed (overlapping repeats are
  merged etc), alternative output (.ace or .gff) can be created and an
  overview .tbl file will be created. Classes that will be displayed in
  the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything
  with 'RNA' in it, 'Simple_repeat', and 'Other' or 'Unknown' (the
  latter defaults when class is missing). Subclasses are plentiful. They
  are not all tabulated in the .tbl file or necessarily spelled
  identically as in the repeat files, so check the RepeatMasker.embl
  file for names that can be parsed into the .tbl file.

Dear, those numbers in the annotation tables not accurate, it's totally different.  i produced these tables to show an example of what I'm looking for, my repeat roughly have all type of repeats in my library file . and all the headers are suits repeatmasker analysis, ready for annotation but when u use the -lib flag the tool assume the analysis is homo and the .tbl is limited, i would like to run the annotation and obtain the second table

I never used repeatmodeler

The approach you have mentioned is beneficial if you want to combine the new de-nove library obtained with those already discovered and deposited earlier in repbase, thank you. But my question here to get more annotation details in.tbl file

I get it now, you asked about the labels, not about the numbers. I added info about how RM annotates them

Dear thank you for the update, but I have edited all my headers earlier corresponding to the order of repeatmasker library documents, but the table still assume to be homo and limited

RepeatMasker comes with a buildSummary.pl script.
You can input the .out file from your RepeatMasker run and there is a flag to enter your species
e.g.
buildSummary.pl -species Vertebrates rmrun.out > mysummary.tbl

In the mysummary.tbl, the summary will not assume your run to be Homo (default).


text:
select unique names while removing the duplicates from a column
I have a file with Ensembl gene and Transcript IDs :

I want to retain only the first entires from both the column so that the output looks like this:
ENSG00000003137 ENST00000001146
ENSG00000003402 ENST00000309955

and ignore the rest.
Have tried awk awk '!a[$1$2]++' and versions of it. Nothing works.
Kindly suggest some way to get the desired output.
regards.

Keep it simple:
awk '!a[$1]++' file.txt

If you have a different command line separator you need to specify with the -F option.
Space
awk -F' ' '!a[$1]++' text.txt

Tab
awk -F'\t' '!a[$1]++' text.txt

If you don't have the first column in order maybe should be a good idea sort first.
sort text.txt | awk -F' ' '!a[$1]++'

The name of the array in this case is a. There is a better explanation in this post.

That sound like a job for a perl one-liner:
perl -ane 'unless(exists $seen{$F[0]}){print $_; $seen{$F[0]}++}'

Explanation:
The input is read line-wise and printed back to the command line if the ID in the first column has not been seen before.

perl -ane ' reads input line-wise and splits fields to array @F
unless(exists $seen{$F[0]}){ checks if ID exists in hash %seen and continues if that is not the case
print $_; prints the original input to STDOUT
$seen{$F[0]}++}' adds the ID to the hash %seen


text:
Scientific names of taxids with kraken 2.0.8 beta, command kraken2 --use-names error
I am running kraken 2.0.8 beta to do taxonomic classification of micorbiota data coming from a MinION run based on the 16S gene. To do so I am using the pre-bult kraken SILVA  database. I managed to get the standard kraken output format but now I would like to replace the taxonomic ID column with  scientific names. According to the manual (https://ccb.jhu.edu/software/kraken2/index.shtml?t=manual) I should run the command:
kraken2 --use-names
but it gives the error:
command kraken2 not found
I also tried variants such as:
./kraken 2 --use-names
but still it gives an error.
I also tried without kraken2 and thus,
--use-names
but it gives this error:
--use-names command not found
Do you have any idea how I can proceed?
In the picture you can see that everything works to obtain the standard kraken output for barcode 76. But if I try to implement the --use-names command then it does not work

I solved the problem: it was a matter of the position of the --use-names command. By changing position (see picture) it worked.

What is the output when you run which kraken2 in your command line?

Simply nothing happens

Install kraken2 and add the executable to your path as per https://ccb.jhu.edu/software/kraken2/index.shtml?t=manual#installation , or use a complete absolute path to your kraken binary when running it e.g.: /path/to/your/kraken/install/bin/kraken2 --use-names <rest of your command>

I think that I am already using the absolute path. you can see it in the picture I added to my post


text:
Get start and end coordinates per chromosome
I will use GATK for SNP calling (HaplotypeCaller). I need to feed the interval file in the command, otherwise I get errors (even though I want to use the whole genome, not a subset  - btw it's not the point of the question, so will not elaborate).
The interval file for GATK can be for example a bed file, with 3 columns: chr_name  chr_start  chr_end. I don't have this file, but have the genome and the reads.
To obtain the intervals per chr/scaf, I proceeded this way:

converted the alignment to interval file, using bam2bed    (BEDOPS)
extract the columns $1,2,3 (using cut)
get the min & max coordinates per chr/scaff with a small R script:

bd <- read.table('all.bed', h=FALSE)
bd_min <- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max <- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord <- merge(bd_min, bd_max, by='bd$V1')
write.table(bd_coord, file='coords.bed', sep='\t', quote=FALSE, col.names=FALSE, row.names=FALSE)

These should thus represent the mapped intervals, conceptually same as not feeding the intervals at all (= considering all the sequences).
The ideas for this:

in the bed file: $1 is be the feature, $2 the start position, $3 the end position
so, I need to get the min of $2 for each unique item in $1(for the start coord) and the max of $3 for each item in $1 (for the end coord)

The output looks plausible:
$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560

My question here:

is this a correct way to proceed?
is there a standardized way to perform this?

My main concern is that I will have to re-build this "coordinate" files for every GATK run, because I think the coordinates per chr/scaf will shift, even slightly, at each GATK run with different datasets.
My other minor concern, is that this approach is a bit slow (given the 160M lines in the current bed file), so an unix-tools solution will also be accepted. I tried to compute myself the max/min using awk (e.g. awk '$3>max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'), but I get different results than with the R approach...

You can try this, since you have the bam file, all of them are aligned to the same reference genome (I hope)(and apologies for the awful awk.. been a while since I used it):
samtools view -H some.bam | awk '{if($1~/@SQ/){print $2"\t1\t"$3}}' | sed 's/[SL]N://g' > coords.bed

For the human genome I got this:
chr1    1   195471971
chr2    1   182113224
chr3    1   160039680
chr4    1   156508116
chr5    1   151834684
chr6    1   149736546
chr7    1   145441459
chr8    1   129401213
chr9    1   124595110

What I did is to get the chromosome length, and make every line start at 1.
If you use R, you can do:
library(Rsamtools)
CL = scanBamHeader(<bam file>)[[1]]$targets
write.table(data.frame(names(CL),start=1,end=CL),
"coords.bed",row.names=FALSE,quote=F,col.names=FALSE,sep="\t")

thanks for your input @StupidWolf. yes I'm using the same reference, to which i also aligned all the reads. it helps the question, on how to get chromosome coordiantes - although not solving the GATK error,  still getting "Input files reference and reads have incompatible contigs: No overlapping contigs found.")

Weird... are you providing the fasta file with the same header? grep ">" <your.fasta> and samtools view -H  | grep "@SQ" ... you should see the same chromosome names

yes, weird, agreed. and yes, i verified that they contain the same fasta IDs, which I also see in the bam header... I opened an issue on GATK forum, as self-comment to my question: https://gatk.broadinstitute.org/hc/en-us/community/posts/360057671431-GATK-HaplotypeCaller-errors-with-or-without-the-L-interval-file-

I upvoted your answer since it's what I was asking for in the original question. I'll deal with the GATK issue on my own. :up:

Thanks dude.. Can you make a small bam file and provide the fasta via ftp or something? Most likely the GATK people might ask for it. I cannot use GATK right now, but I can try with bcftools..

Not sure if I completely understood your question, but if you want simply an interval of your entire genome  you could do instead:
samtools faidx MyGenome.fasta

Obviously that means that you need access to that, not sure if you have.
The resulting fai file contains a format similar to that:
000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

Where the first column is your chromosome and the second one your length.
I guess you could then just do a awk '{OFS=FS="\t"}{print $1,"0",$2-1}' to get your file.

yes, I do have the genome, and have indexed it. as reported in the question itself, I could obtain (what I understood being) an interval file (just: chrID, start, end). and with @StupidWolf's answer I just tried to expand this interval to start from pos1 for each chr/scaf feed it in GATK. but I still get the same GATK errors (full link in my self-comment to my question) I got for the other interval files. weird, I agree. 100% sure I used one and only ref to map the reads to. I only have that one.

I think the reported error from GATK is clear, the annotation which should look as following:
Bla_scaf1190  830  763
Is wrong as the start is higher then the stop. The question is why you would make an interval file which covers only mapped regions. Isnt the point of intervals to reduce the mapped space ? Otherwise speed should be identical for full space and your intervals

that's one test I tried, and that one error is clear. in the two other options (1.no intervals at all, technically taking the whole genome;  2.intervals for mapped regions of the chr; and I could add what I just tried, so 3.intervals for the whole chr length), I have the other error still persisting "Input files reference and reads have incompatible contigs: No overlapping contigs found"

This probably comes too late, but for the record GATK supports multiple kinds of interval file format. BED is the most common, and there's the Picard format (.interval_list) which is BED-like with a few tweaks. Both require coordinates. However, GATK also supports a simpler format (.list or .intervals) that allows you to just give the contig names (eg chr1) without any coordinates. Sounds like that would satisfy your need. Incidentally, GATK also supports using VCF files as interval file if you want to re-genotype specific regions based on variants located there (ideally paired with the interval padding arguments).
See the GATK doc article for more details:
https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists
(The doc has a few formatting issues but I just notified the team so they'll fix it)


text:
How to generate ultrametric phylogenetic trees?
This question has also been asked on Biostars
I am going to use ultrametric tree for downstream analysis. I would like to see use the ultrametric tree for Ranger-DTL for infer gene family evolution by speciation, gene duplication, horizontal gene transfer, and gene loss.
I have generated rooted species tree with 1000 bootstraps by RAxML. Further, I need to generate ultrametric tree for the same. Therefore, Please suggest me how to generate ultrametric tree from the rooted species tree.
Thank you.

An ultrametric tree is simply a global molecular clock and can be easily generate in PAML using a single evolutionary rate. You take you RAxML tree in phylip format (PAML might do nexus format know), and click the "clock function". PAML will do the rest.
You can alternatively generate it de novo via UPGMA neighbor-joining tree which is found in the MEGAX package (which is probably easier and considered a technical definition of an ultrametric tree).
A global molecular clock (as opposed to a local molecular clock) is quite a specific request, but if you performed a likelihood ratio test, it would identify "clock-like" behaviour.


text:
What are the meanings of these transcript ids?
I have three types of transcripts for Rosa chinensis "Old Blush" homozygous genome v2.0 from GDR genome browser.
RcHm_v2.0_CPg0501991
RcHm_v2.0_Chr0c12g0499591  or RcHm_v2.0_Chr1g0313011
RcHm_v2.0_MTg0498341

Does anyone know what this means. At least "CPg" and "Mtg". I know "Chr" means chromosome. And there is "1" at each and every transcript. But it is not separated using anything, like people do to mention the transcript number. I will be grateful if someone can help.
Thank you.


text:
How to create a .phy file for constructing a phylogenetic tree?
I'm new here. I need to construct a phylogenetic tree from 16s Sequence data (merged 151 bp) and color the branches/nodes based on my experimental conditions. From this blog, it looks like I need a .phy file to make the tree. I want to use Biopython for this since I'm reasonably comfortable with python.
I have a csv file with the 16s sequences, the phylogenetic information, and the "name" of the sequence, and the experimental conditions. How can I make a .phy file from this?
I would appreciate any help.
Thank you.

You simply use Biopython's SeqIO and AlignIO.
The documentation is https://biopython.org/wiki/AlignIO
etc ...
The format you seek is phylip

Thank you for your response. I made a .phy file similar to the one here link (except that in mine, the length of the sequences is 253). I ensured that they all have the same length.  But when I ran  al=AlignIO.read(open("for_tree_top10.phy"),"phylip"), I got "ValueError: Sequences must all be the same length". I'm not sure why, because I explicitly checked that the sequences have the same length. Not sure what to do now :( Any thoughts?

Possibly it is not a programming issue but formating. There is no simple answer to trouble shoot this. Could you ask this as a separate question please? Obviously it is solvable but there's several possibilities and the Moderators will intervene if there are multiple questions within a thread

I have asked the question here link

thanks answered it

Python works great for tying everything together and file manipulations.  I find stepping out of Python, using system calls, best for alignment and tree building.  For example:
os.system("mafft --auto {} > {}" .format("unaligned_16s_seq_in_file.fasta", "alignment_out_file.fasta"))

os.system("raxml -f a -s {} -p 12345 -x 12345 -# 100 -m GTRCAT -n {}" .format("alignment_out_file.fasta", "tree_out_file.tre"))


text:
How can i get data set of illness and not illness people?
I'm reading a virology paper where the authors made a classification for the patients as to whether they had a specific disease outcome or not against the underlying sequence ...

The research focused on two different class of DNA which are Hepatitis B Virus
(HBV) Genotype B nucleotide sequence and Non-HBV nucleotide sequences`

How can I get the data set or how can search to find it ?

How can data sets of sequences date classify between illness and asymptomatic patients?

Please edit your question and add a link to the paper you're referring to. If the paper does not point to the dataset it analyzes, the best way to access the data (if it's not restricted access) is to email the corresponding author. Also, what do you mean by "sequences of illness and not illness people"?

Also, many researchers deposit their data at NCBI GEO: https://www.ncbi.nlm.nih.gov/geo/
...loads of stuff there.

The data set will be available on Genbank, you simply search for the "authors surname and HBV" and their entire data set should be revealed.
The HBV genome is 5Kb, the smallest viral genome so its not alot of data.
I am confused by "Genotype B nucleotide sequence and non-HBV nucleotide sequences". What I suspect is the investigators looked at HBV genotype B in comparison to HBV genotypes A, C to I. Essentially they wanted to assess the relationship between genetic variation and pathogenesis, which is an attempt to assess virus-induced cancer. Obviously I don't know the paper but I'm aware of the basic research front in this area. HBV has a very wierd evolution and I personally wouldn't want to investigate it.


text:
bwt suffix array construction
I found a very good link to generate bwt suffix array in BW based tools.
@user172818 provides a nice code. could anyone tell me how to store SA elements' original position in the original string?


text:
SwissProt Proteins into my Jbrowser
How is it possible to get the SwissProt Proteins into my Jbrowser in the same way like shown below?

Thank you in advance,


text:
Snakemake: Migrating from deprecated cluster.json to new profiles.yaml
I am an avid user of Snakemake. Recently we have been refreshing our pipelines and I saw that a cluster.json file is no longer the recommended way to store the cluster configuration.
I used to start my pipelines like this:
snakemake --cluster-config cluster.json --cluster "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}" --jobs 245 --use-conda --rerun-incomplete --keep-going --latency-wait 60

With the cluster.json looking like this:
{
    "__default__" :
    {
        "cluster"   : "qsub",
        "jobs"      : 70 ,
        "nodes"     : 1,
        "ppn"       : 1
    },
    "rsem_quantify" :
    {
        "ppn"       : 8
    },
    "star_map_human" :
    {
        "ppn"       : 4
    },
}

Now, I understand that some of this goes in to the new profiles.yaml like this:
cluster: qsub
jobs: 70

But how do I set the ppn values per rule? And how do I set the other configuration options into profile.yaml (like: snakemake "qsub -W umask=002 -d cluster_outputs  -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {rule}", etc.)
Any examples are appreciated, https://github.com/snakemake-profiles is not really helping

Just wanted to note that I asked the same question in stackoverflow. No  good answer yet though.

It looks like you can just bundle your cluster config in its own YAML file, and reference that from within a profile YAML.  Here's a clumsy little example using a fake job submission script: https://github.com/ressy/example-snakemake-profiles

The current documentation makes it sound like cluster configuration support is going away and being replaced by profiles:

While still being possible, cluster configuration has been deprecated by the introduction of Profiles.

...but from those example profiles, that might be misleading; here's an example template for a profile YAML that references a separate cluster-config YAML:
https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/config.yaml
...
cluster-config: "$((INSTALDIR))/overwrite_cluster_param.yaml" #abs path
...

That cluster-config YAML has the same default and per-rule entries just your JSON version does:
https://github.com/Snakemake-Profiles/generic/blob/master/%7B%7Bcookiecutter.profile_name%7D%7D/overwrite_cluster_param.yaml
...
# specific parameters for certain rules, which need more time/memory

#run_assembler:
#  queue: bigmem
#  time: 100
#   threads and memory definen in config file

So I think you can just set up a pair of YAML files accordingly including options for your rsem_quantify and star_map_human rules.  From my own test it looks like it handles the cluster argument templating just as it did before when using a profile.  (If this is correct we should nag them to make the new documentation more explicit though.)


text:
Is there a "definitive" database for Mobile Genetic Elements?
This isn't quite my speciality, but I've been tasked with finding any evidence of Mobile Genetic Elements (like transposable elements, retrotransposons, bacteriophages, etc.) within some WGS samples.
Is there a "definitive" database used with sequences? How "species-specific" would this be?
Furthermore, are there standard algorithms used to find these MGEs? I suspect we don't have a comprehensive database of all possible MGEs, but it may be possible to query against the sequences we do indeed have (or detect any repetitive elements within the genome).

Would advise doing Denovo identification method structure-based. Then use these public databases (  RepBase, Dfam ) to look for high similarity or homology exploration, this possibility if your species, the repeats not yet discovered, or deposited in the databases.
If the species well-investigated earlier and the mobile element recognised, I assume it will be a straightforward task using species-specific search via Repeatmasker tool with the -species flag to obtain the MobileElement annotation in the genome.
Let me know if you need any data from Repbase related to your species, i have an entree.

Thanks! This is a great answer. For how long has RepBase been subscriber-only?

I think the middle of last year

Back in the Day™ RepeatMasker[1] was the go-to tool everyone used for TE annotation in eukaryotes. At the time, the RepBase database from GIRI[2] (license required) was the de facto database recommended for RepeatMasker users. I know de facto isn't necessarily definitive, but it may have been the only public database supported for some time.
It looks like more recent efforts have focused on developing the Dfam database, mirrored after the design of databases like Pfam and Rfam.

[1]Smit AFA, Hubley R, Green P (2013-2015) RepeatMasker Open-4.0, http://www.repeatmasker.org.
[2]Jurka J, Kapitonov VV, Pavlicek A, Klonowski P, Kohany O, Walichiewicz J (2005) Repbase Update, a database of eukaryotic repetitive elements. Cytogentic and Genome Research 110:462-467, doi:10.1159/000084979.

Thanks for the answer. So, if I understand correctly, the strengths of the Dfam database over RepeatMasker is that it's more extensive?


text:
How I change a function not to need an input
There is a function does need an input as selected genomic site across genomic regions. But my genomic region bed file already is a subset of interesting site (like promoters). I want this function to work not complaining about genomic site
However when I leave genomic site says
> prepare_elements_regions_tf <- prepare_elements("regions.bed",
+                                                 window_size = 50000)
Successfully read elements
 Error in utils::read.delim(sites_file, stringsAsFactors = FALSE, header = FALSE) :
  argument "sites_file" is missing, with no default

The whole function is here
https://github.com/wenmm/ActiveDriverWGS/blob/master/R/ActiveDriverWGS.R
But part of function that nags is this
prepare_elements = function(elements_file, sites_file, window_size = 50000, mc.cores = 1, recovery_dir = NA) {

    try_error = "try-error"

    element_coords_filename = paste0(recovery_dir, "/ADWGS_elements_coords_recovery_file.rsav")

    element_coords_1 = NULL

    load_result = suppressWarnings(try(load(element_coords_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        element_coords_1 = prepare_element_coords_from_BED(elements_file)

        if (!is.na(recovery_dir)) {

            save(element_coords_1, file = element_coords_filename)

        }

        cat("Successfully read elements\n")

    } else {

        cat("Recovered elements\n")

    }

    element_coords = coords_sanity(element_coords_1)

    element_3n_filename = paste0(recovery_dir, "/ADWGS_element_3n_recovery_file.rsav")

    element_3n = NULL

    load_result = suppressWarnings(try(load(element_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        element_3n = add_trinucs_by_elements(get_3n(element_coords))

        if (!is.na(recovery_dir)) {

            save(element_3n, file = element_3n_filename)

        }

    }

    window_coords_filename = paste0(recovery_dir, "/ADWGS_window_coords_recovery_file.rsav")

    window_coords = NULL

    load_result = suppressWarnings(try(load(window_coords_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        window_coords = coords_sanity(get_window_coords(element_coords, window_size))

        if (!is.na(recovery_dir)) {

            save(window_coords, file = window_coords_filename)

        }

    }

    window_3n_filename = paste0(recovery_dir, "/ADWGS_window_3n_recovery_file.rsav")

    window_3n = NULL

    load_result = suppressWarnings(try(load(window_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        window_3n = add_trinucs_by_elements(get_3n(window_coords))

        if (!is.na(recovery_dir)) {

            save(window_3n, file = window_3n_filename)

        }

    }

    raw_site_coords_filename = paste0(recovery_dir, "/ADWGS_raw_site_coords_recovery_file.rsav")

    site_coords_filename = paste0(recovery_dir, "/ADWGS_site_coords_recovery_file.rsav")

    raw_site_coords = NULL

    site_coords = NULL

    sites_load_result = suppressWarnings(try(load(site_coords_filename), silent = TRUE))

    if (class(sites_load_result) == try_error) {

        load_result = suppressWarnings(try(load(raw_site_coords_filename), silent = TRUE))

        if (class(load_result) == try_error) {

            raw_site_coords = utils::read.delim(sites_file, stringsAsFactors = FALSE, header = FALSE)

            colnames(raw_site_coords) = c("chr", "starts", "ends", "site_info")

            raw_site_coords$chr = gsub("chr", "", raw_site_coords$chr, ignore.case = TRUE)

            raw_site_coords$chr = paste0("chr", raw_site_coords$chr)

            if (!is.na(recovery_dir)) {

                save(raw_site_coords, file = raw_site_coords_filename)

            }

            cat("Successfully read active sites\n")

        } else {

            cat("Recovered active sites\n")

        }

        if (!is.null(raw_site_coords[[1]])) {

            # first only keep sites that overlap with element coordinates

            site_gr = GenomicRanges::GRanges(raw_site_coords$chr,

                IRanges::IRanges(start=raw_site_coords$starts, end=raw_site_coords$ends))

            element_gr = GenomicRanges::GRanges(element_coords$chr,

                IRanges::IRanges(start=element_coords$starts, end=element_coords$ends))

            rm(raw_site_coords)

            overlaps_found = GenomicRanges::findOverlaps(element_gr, site_gr)

            site_here_gr = site_gr[S4Vectors::subjectHits(overlaps_found)]

            rm(site_gr, element_gr, overlaps_found)

            site_coords = do.call(rbind, parallel::mclapply(1:nrow(element_coords),

                get_sites_per_element, element_coords, site_here_gr, mc.cores=mc.cores))

            rm(site_here_gr)

            site_coords = coords_sanity(site_coords)

        }

        if (!is.na(recovery_dir)) {

            save(site_coords, file = site_coords_filename)

        }

    } else {

        cat("Recovered active sites\n")

    }

    rm(sites_load_result, element_coords)

    site_3n_filename = paste0(recovery_dir, "/ADWGS_site_3n_recovery_file.rsav")

    site_3n = NULL

    load_result = suppressWarnings(try(load(site_3n_filename), silent = TRUE))

    if (class(load_result) == try_error) {

        site_3n = get_3n(site_coords)

        if (!is.na(recovery_dir)) {

            save(site_3n, file = site_3n_filename)

        }

    }

    list(element_coords=element_coords_1, element_3n=element_3n, window_coords=window_coords, window_3n=window_3n, site_coords=site_coords, site_3n=site_3n)

}

How I get rid of nagging part?

What have you tried?

I tried to edit the function by elimination any line goes with site fine but function is not working at the end

Think about it this way - your elements_file already has only the sites of interest and hence does not need to be subject to a subset operation. In other words, if you were to have a sites_file, it would have the same content as the elements_file.
Can you see where I am going with this?

I think you are stating to edit the function for any site file requisition

You don't need to edit the function. sites_file would have the same content as the elements_file. What does this mean in terms of passing file names to the parameters? Are they different files?

sites_file is a BED4 separate input file but as the content already in element_file, I don't provide that. Element_file is a BED12

So, sites_file should contain fewer columns and possibly (but not necessarily) a subset of the regions in elements_file, correct? If they contain the same regions, is it not straighforward to get the BED4 content from a BED12 file? You literally need 4 columns of an existing file for a new file.

You right contains 4 columns for chromosome start and a more column that I don't know what could be like this  chr17 7565257 7565267 site1. My problem is I don't site column is which part of BED12 file with 12 columns

Please invest some effort in learning the columns in various BED formats. I literally googled "BED4" and got this: https://bedtools.readthedocs.io/en/latest/content/general-usage.html. This should tell you how to get to BED4 from BED12.


text:
Can I get longer contigs by changing MEGAHIT settings?
My current settings are megahit -r input.fq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output. I picked the 'presets meta-large' because I am dealing with a quite large and diverse metagenome out of which I only need some genes.
I've got a contig of interest ~19kb, with a depth of 16. Could I make this contig larger with more relaxed settings? Help would be appreciated!

You can use the option --k-step 10 and --kmin-1pass for low-coverage datasets.


text:
How to filter a SAM file by a bed file?
I have a large BAM file and I want to filter it by just what's defined in a bed file, and write it to a new file.
Is it possible to do this with samtools, and if so, how? Is there a better tool I could use to do this?

Did you look at the options for samtools view?

´samtools viewhas filtering options for BED files as reference. Typesamtools view` and check available flags.

Hi, i.think your question needs more specific details

Thanks to the helpful comments I figured it out:
samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam

-b = Output as bam file
-h = include header
-L = only output alignments overlapping the input BED file

You can also do this with bedtools intersect:
bedtools intersect -abam input.bam -b bedfile.bed > output.bam

This should be identical to the samtools view answer. One further feature though is you can output all reads that don't overlap with the regions in bedfile.bed by adding the -v flag.


text:
What's the difference between Procheck and Whatcheck?
I'm currently using Procheck* to check a protein backbone. The process is:

Determine backbone phi, psi, and omega angles
Compare to known angle distributions to identify outliers, and angles in disallowed regions

Is Whatcheck** equivalent for this task?
*, checks the stereochemical quality of a protein structure, producing a number of PostScript plots analysing its overall and residue-by-residue geometry. It includes PROCHECK-NMR for checking the quality of structures solved by NMR.
**, Derived from a subset of protein verification tools from the WHATIF program (Vriend, 1990), this does extensive checking of many sterochemical parameters of the residues in the model.


text:
STAR quantMode vs RSEM vs Kallisto
I recently discovered this Snakemake pipeline for RNASeq that uses STAR's quantMode to quantify gene expression for DESeq2 differential expression analysis.
In the past I've always seen the workflow as:

STAR alignment
RSEM/HTSeqCount quantification
DESeq2 analysis

More recently I've also seen:

Kallisto/Salmon pseudo-alignment quantification
DESeq2/Sleuth analysis

I tried searching online for differences between STAR quantMode vs other quantification algorithms but couldn't find many details. What are the benefits and drawbacks of using STAR quantMode vs RSEM/Kallisto/Salmon?

STAR quantMode (GeneCounts) essentially provides the same output as HTSeq-Count would, ie. number of reads that cover a given gene. This is the most simple measure of expression you could get from RNA-seq data. Kallisto and Salmon utilize pseudo-alignment to determine expression measures of transcripts (as opposed to genes). RSEM uses some algorithm to determine isoform fractions, but still based on alignments. (Note that I don't have any personal experience using RSEM.)
It ultimately depends on the research question which method will be most appropriate for you.
read counts
pros:

relatively easy to understand concept-wise
can be used to generate many different expression measures (FPKMs, TPM, CPM, etc.) and is, thus, more versatile in down-stream analyses

cons:

you can't distinguish isoforms (from the counts themselves)
requires alignment and, thus, more time to obtain

salmon/kallisto:
pros:

faster (due to lack of alignment step)
quantifies (known) isoforms

cons:

expression values are more abstract, less easy to explain
gene-level analyses will require some extra work

RSEM
Likely lies somewhere in between, still requiring alignments will make it slower. But it does give you expression values for isoforms.

A shorter version...STAR's is convenient, in that it gets done right alongside the alignment, but RSEM and Kallisto are smarter about dealing with ambiguous reads.  So I'd say your results from either of those are superior.
Personally, I use STAR for alignment, then give the bam to RSEM.  STAR's gene count is good for a reality check to make sure that I got the strandedness right, and RSEM's gene level expected counts are good for DESeq once they are rounded.

You can also feed the transcriptome coordinate BAM file generated by STAR directly to salmon to quantify abundances if you want to rely on the STAR alignments.  As an added bonus, if you process the STAR alignments with salmon, you need not disallow indels (gaps) or soft clips in the alignments (--quantTranscriptomeBan Singleend instead of IndelSoftclipSingleend).  This can be useful if your transcriptome has variation with respect to the reference.


text:
Protein ligand docking: how to convert <protein>.pdb to <protein>.maps.fld?
Hello I'm helping to develop a cloud docking tool for screening compounds, similar to Swissdock but with mass throughput and GPU optimizations.
Specifically helping screen existing drugs against coronavirus proteins (couple structures solved, rest TASSER).
I'm planning to use GPU optimized AutoDock-GPU , which takes in .maps.fld

I know you can use autogrid to select the bounding box and generate the .maps.fld, but I've been unable to figure out the workflow.
for preliminary screening I want to search the whole protein without specifying a bounding box.
Is there a script for converting protein.pdb to .maps.fld?

Collaborators welcome.
Thanks!

Hi @PaulShen please do remember to mark the answer 'accepted'. This really does help the site. The answer is very thorough

There are a few steps in between pdb and .maps.fld.
Here is the list of several scripts that can do tasks for you that you downloaded with autodock MGLTools: http://autodock.scripps.edu/faqs-help/faq/where-can-i-find-the-python-scripts-for-preparing-and-analysing-autodock-dockings. Look at the prepare_ files. Also note that the scripts come with a pre-bundled python, but you can install these scripts with conda.
Here are the basic order of things:
prepare_receptor4.py -r protein.pdb
prepare_ligand4.py -l ligand.mol2
prepare_gpf4.py -l ligand.pdbqt -r protein.pdbqt -y
autogrid4 -p protein.gpf
prepare_dpf4.py -l ligand.pdbqt -r protein.pdbqt

partial charge
Your PDB needs partial charges so first you convert it to PDBQT file.
This is true for your protein and your ligand. That is no complicated parameterisation for your ligand —which is not necessarily a good thing!
In this step it is essential to correct protonation if absent each have their own flags.
grid
The you make a grid file (gpf), which contains your grid. To make a box that spans the whole protein simply use the unit cell dimensions form the PDB. For an mmCIF, if you open the header dictionary you need _pdbx_struct_oper_list.matrix using MMCIF2Dict from biopython's PDB submodule. For a PDB file, just search for the CRYST1 line (cf. format).
Once you have have the grid you can run autogrid. This step will create the map.fld
Why a big box is bad for a VS
However, all those who take this approach opt for a box that is huge: as there is no solvent having a huge bounding has little penalty. However, this is considered a very bad strategy for a virtual screen. This not because it uses up computer resources as these are dirty cheap. Structural knowledge is important and targeting the active site of an enzyme stops it, while binding to the surface does nothing unless it's an interface. So manually reading what the stuff are and choosing the box wisely will save you a lot of time downstream in analysis. There is after all the saying "one week in the lab will save you an hour in the library"...
Protein and ligand
Lastly, merge the protein with the small molecule in a dpf file.
Caveat against docking with models
It is not a good idea dock against models. Swissmodel (unless virtually identical), I-Tasser, Phyre or EVFold models. Docking is very sensitive to small structural difference that may have resulted from the modelling, so docking to models is highly discouraged. For Coronavirus protein, dock against SARS, which is highly similar. Or at the very least pay strong attention to the I-Tasser C-scores and discard the bottom 2/3 of their models.
Note about coronavirus
Your objective of docking coronavirus may be a bit too late however.
I did a wee interactive summary of the literature about the solved structure which is the protease. Proteases are very easy to rationally design drugs for.
There are already:

The Covid Moonshot project is an open collaboration between a fragment-screening X-ray facility (Diamond XChem) and many researchers from various disciplines (myself included).
a paper that has an empirically validated coronavirus specific ligand
many many manuscripts about virtual screens in bioarxiv/chemarxiv.
licensed HIV protease inhibitors lopinavir and ritonavir that were initially suspected to be effective —a clinical trial

Also you may want to check out Galaxy project, a tool to run pipelines for genomics, but have started moving into biochemistry from genomics. Eg. https://covid19.galaxyproject.org/cheminformatics/.

Thanks Matteo for your detailed explanation and literature references! It's great to see promising candidates in the pipeline. Our aim is mainly on the tool--a modern GPU optimized Swissdock-like web service for quickly screening thousands of existing drugs. Then automatically suggesting candidates not only on the affinity but also clinical plasma concentration, bioavailability and side effects. The results wouldn't be novel for coronavirus though good if we reproduce existing literature.

A wee suggestion. I understand that most users may not be familiar with Lipinski's rule of five or normal predicted chemistry metrics, such as logP, but these are likely more useful than pharmacology metrics, especially as the latter are very rarely available. Plus the Crippen contribution to logP is always found intriguing to non-compchemists, eg. Bioinf SE question about lipophilicity

That makes sense. Thank you for the advice.

The Scripps link is broken the following provides an avenue to find it https://autodock.scripps.edu/resources/


text:
Prevent Snakemake from changing wildcards: sample for every rule
EDIT: I posted this as a bug, perhaps wait for the result of that? github
All,
I am using Snakemake to build a pipeline. To define input and output names, I constructed them based on a filename base (from the samples table), a read indicator (like R1, R2) and the fastq file extension, like so:
samples = pd.read_csv('samples.tsv', sep='\t').set_index('samples', drop=False)
reads = ['R1', 'R2']
fastq_file_extension = _001.fastq.gz

rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

However, during pipeline execution already I can see strange things happening, for this rule it reports:
[Fri Feb 21 12:39:57 2020]
rule fastqc:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz
output: qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.html, qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.zip
jobid: 10
wildcards: sample=0053_P2017BB3S20R_S2_R, read=1

As you can see, Snakemake pulled the R part of the read to the sample and set the read to 1 (should be R1). But this still works because it produces output files as Snakemake expects.
But the next rule, trim_galore_pe:
rule trim_galore_pe:
    input:
        [os.path.join(config['Directories']['fastq_raw_dir'], '{sample}' + read + fastq_file_extension) for read in reads]
    output:
        fastq_read1 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.1_val_1.fq.gz'),
        fastq_read2 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.2_val_2.fq.gz'),
        report_read1 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.1.fastq.gz_trimming_report.txt'),
        report_read2 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.2.fastq.gz_trimming_report.txt')
    params:
        extra="--retain_unpaired"
    wrapper:
        "0.49.0/bio/trim_galore/pe"

This goes wrong, during pipeline execution Snakemake reports:
[Fri Feb 21 12:40:37 2020]
Error in rule trim_galore_pe:
jobid: 7
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
conda-env: /home/nlv24077/projects/snaqs/snakemake_envs/1b94430b
cluster_jobid: 11904662.osios
Error executing rule trim_galore_pe on cluster (jobid: 7, external: 11904662.osios, jobscript: /home/nlv24077/temp/RNASeqPairedEnd/.snakemake/tmp.ooypuklc/snakejob.trim_galore_pe.7.sh).
For error details see the cluster log and the log files of the involved rule(s).

This is understandeble because the file fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz does not exist, it should be fastq_trimmed/0053_P2017BB3S20R_S2.1_val_1.fq.gz, somehow Snakemake added an extra underscore, which you can see in the execution output, under "wildcards" (0053_P2017BB3S20R_S2_).
[Fri Feb 21 12:39:57 2020]
rule trim_galore_pe:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz, fastq_raw/0053_P2017BB3S20R_S2_R2_001.fastq.gz
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
jobid: 7
wildcards: sample=0053_P2017BB3S20R_S2_

So how can I prevent Snakemake from changing the samples and read wildcard in seemingly arbitrary ways?

Placing two wildcards directly behind each other will make it impossible for snakemake to determine where one wildcard ends and the other begins. So it is simply taking the longest possible match for the first wildcard (sample) and then one character for the second (read). This will also cause the _ to be included in the wildcard value sample.
I'm assuming that this same pattern is used for the inputs of another rule to refer to the trim_galore_pe output. Hence why the _ is included in the wildcards there.
You might be able to fix this, simply by adding a _ in-between {sample} and {readgroup} in your first rule (and any other rules using this same pattern):
rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}_{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

I did this but Snakemake conjures up a new _ out of nowhere:
Building DAG of jobs...
MissingInputException in line 2 of /rst1/temp//RNASeqPairedEnd/rules/trim_galore_pe.smk:
Missing input files for rule trim_galore_pe:
fastq_raw/0053_P2017BB3S19R_S1__R1_001.fastq.gz
fastq_raw/0053_P2017BB3S19R_S1__R2_001.fastq.gz

But why is is impossible for Snakemake to separate them? They are not in {} for nothing right? I don't understand why Snakemake would want to change my wildcards from anything it gets from the samples table. What would be the use? It seems to me that this is just a bug? Or am I really fundamentally misunderstanding Snakemake?

Regarding my first comment, one might think that Snakemake now added and extra _ so I can remove the existing one perhaps? But no, it then removes both:
Building DAG of jobs...
MissingInputException in line 2 of /rst1/temp/RNASeqPairedEnd/rules/trim_galore_pe.smk:
Missing input files for rule trim_galore_pe:
fastq_raw/0053_P2017BB3S19R_S1R2_001.fastq.gz
fastq_raw/0053_P2017BB3S19R_S1R1_001.fastq.gz
So it's either 0 or 2, not 1, as one would expect.

So I guess I am misunderstanding this completely, When I change the input of trim_galore_pe, it complains of missing files, but when I remove the input entirely it complains that the fastqc rule is missing input files, even though I did not touch that rule. Is there any way to more clearly see where errors originate exactly? I feel completely lost now.

@Freek Could you share the entire Snakefile? It's quite hard to see what's going on exactly without it.

Ok, I think I have found the reason, I was using code as well to create the samples_table from the same config and it was adding and removing the _ there as well. So this method does work, my apologies. It remains unhandy that I can't directly use wildcards together, But I could perhaps define a separator in the config file to solve this and keep flexibility..

I think you want to prevent snakemake from interpreting your wildcard strings as regular expressions. For this, escape each wildcard value using wildcard_constraints. E.g., towards the top of your Snakefile add:
wildcard_constraints:
    read= '|'.join([re.escape(x) for x in reads]),
    sample= '|'.join([re.escape(x) for x in samples])

(Personally, I use the wildcard_constraints directive as above quite liberally to constrain every wildcard)

Does this mess with how you can use rules.fastqc.output.html as the input for other rules? It gives me a lot of errors for a downstream step (multiqc).

As far as I know, this doesn't cause problems. If you don't show the error it's difficult to tell what's going on.