text:
Create GFF3 feature exporter
Problem Statement

The GFF3 format is a commonly-used one in bioinformatics for
representing sequence  annotation. You can find the specification
here:
http://www.sequenceontology.org/gff3.shtml
The genome and annotation for Saccharomyces cerevisiae S288C is on the
class server here:
/home/jorvis1/Saccharomyces_cerevisiae_S288C.annotation.gff

Note that this same file has both the annotation feature table and the
FASTA sequence for the  molecules referenced. (See the '##FASTA'
directive in the specification.) Within the feature table another
column of note is the 9 th , where we can store any key=value  pairs
relevant to that row's feature such as ID, Ontology_term or Note. Your
task is to write a GFF3 feature exporter. A user should be able to run
your script like this:
$ export_gff3_feature.py --source_gff=/path/to/some.gff3 --type=gene --attribute=ID --value=YAR003W

There are 4 arguments here that correspond to values in the GFF3
columns. In this case, your  script should read the path to a GFF3
file, find any gene (column 3) which has an ID=YAR003W (column 9).
When it finds this, it should use the coordinates for that feature
(columns 4, 5 and 7)  and the FASTA sequence at the end of the
document to return its FASTA sequence. Your script should work
regardless of the parameter values passed, warning the user if no
features were found that matched their query. (It should also check
and warn if more than one  feature matches the query.)
The output should just be printed on STDOUT (no writing to a file is
necessary.) It should have a header which matches their query, like
this:

gene:ID:YAR003W .... sequence here ...

As an extra challenge, you can format the sequence portion of the
FASTA output as 60-characters per line, which follows the standard.

So far I have this below. I'm not even sure if this is correct to start with:
#!/usr/bin/env python3

for line in open("Saccharomyces_cerevisiae_S288C.annotation.gff"):
   line1 = line.rstrip()
   if line.startswith("#"):
       continue
   column = line1.split("\t")
   if len(column) != 9:
       continue
   id = column[8]
   type = column[2]

Can anyone help me with this? I'm totally lost. Thank you!

Was there any restriction about using already developed packages for this purpose? If not : https://biopython.org/wiki/GFF_Parsing , you can wrap your script around with this plus argparse to work from the console: https://docs.python.org/3/library/argparse.html

@Avamys please add this as a separate answer. Comments should only be used for transient information that can help improve or supplement questions or answers. See here for more information.

Your answer looks like a good start. You'll probably want to use sys.argv, argparse, or another argument parsing tool to read in options from the command line, then compare them to the fields (e.g. the "type" argument gets compared with the type column).
Given that attributes can be arranged in an arbitrary order, you'll probably need to read all of them from each line and store in a dict in order to determine whether the target is present. A regular expression search is another way to do that.
There are a few libraries available (e.g. from biopython) for reading and subsetting FASTA records.


text:
seqff script got trouble after samtools treatment
I'm trying to study fetal fraction by using a seqff script from here. The instruction says that I need a headless sam file to work so after alignment with bwa, I process my sam file by using samtools:
samtools view -S -q 20 -F 4 target.sam > treated.sam

This process also remove the @SG header since there is no reference file. I also use grep to remove header of the target.sam got after alignment for comparison. However, when I compare, I noticed that ENet calculation on the treated.sam returned with NA results, so half of the method to calculate fetal fraction of this script was gone. I have tested it in multiple cases and it seems that without using samtools, the script worked well but with samtools Enet calculaton not worked. Can someone tell me what is wrong with the script? I'm using R 4.0 and all modification I made to the script is on the args since it cannot recognize args I input with anymore.
p/s: this is few lines of headless target.sam made by using grep
W2003120038L1C001R0050000009    0   chr10   6490952 0   50M *   0   0   GAGCAGAAAAACTGGAAACTCTAAAAAGCAGAGTGCCTCTCCTCCTCCAA  IBGIBAB<ABBCA9IB@BHDIA=ABBB8IB9A;DGII@I8IIAII<IDCA  XT:A:R  NM:i:0  X0:i:151    XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26

And this is from treated.sam, the sam file be processed by samtools
W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26
W2003120038L1C001R0050000212    0   chr4    157022243   37  50M *   0   0   AAAAATAAAAAAGCATTCTGCCCAGATCTCTCAACATTATGTTATTTCTT  4<AA?D2,AAA@IF@DDIDIIDIA2AD=DBAF@BIADA<DICD@DADFD@  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000258    0   chr13   106093727   37  50M *   0   0   AATGGTTACTTTATGTTTATATTTCAACAGGGAATATTTGAAATATGGCA  97?IIDA*FDAD<@H@AD<@<?AAG8??,I;I@:DADDAG7@@7(;IH6?  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7C42

@gringer I have added few lines of both file on my post. Also I have to perform quality filter to filter out bad alignments results

If the grep one works, and the non-grep one doesn't, my first step would be to exclude any confounding factors (like quality filtering) to try to get them both to look the same (which they don't, at the moment). Is there any improvement if you don't include the quality filter (either in the files looking similar, or the script working better)?

If you think there is a formatting problem, you might need to put small subsets of the samtools file through until you can identify the problematic line.

@swbarnes2 thanks, do you know any command or program that can check for the different? Because the each sam file is more than 4GB, finding what is different by compare with text editor is impossible to say at least

One you figure out the problem in one file, you should be able to fix it for the others.

source: Stack Exchange

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

sort: votes


text:
RNAseq DE comparison with samples of different read length
I have RNA samples but with different read lengths (Eg, HiSeq 2x125 and NovaSeq 2x150bp data). I would like to do DE analysis on these samples. What do you recommend? Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or any other suggestions.

Does this answer your question? Normalize RNA seq data from multiple runs for expression analysis

You don't need to do anything, this is taken care of by the linear model - if you don't have any confounding effects.
If you do have confounding effects, e.g. all control samples on HiSeq and all treated on NovaSeq there's not a lot you can do anyway, other than be honest when you present your results.
Taken care of by the model means adding an extra 2 level factor for sequencer by sample, you need to do this anyway to take care of the batch effect between sequencer runs.

Yes, I would simply trim the longer one to the shorter read length, then map the data as usual. I doubt though that it makes a big difference on the global scale but the processing step is simple and it makes sure you have no mappability bias. If you do not have any other confounding effects then you can do a DE analysis without any extra covariates, I think it is always desirable to keep things as simple as possible and trimming is the simplest solution.

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


text:
Which database is good to store (metagenomic) relative abundance data?
The simple question is: how to store relative abundance data in a database? Which database to choose?
I am talking here about metagenomics data, relative abundance of species or some features. Each feature has single abundance value. For each sample, there are 100's to 1000' features, and the number of unique features between samples keeps growing fast with new samples. Target number of samples ranging 1000-10000.
The most lame way of handling that would be some big array (e.g. sampleID x featureID) but that would not be very convenient, as it will be holding mostly zeros (for features not found in given sample).
I wanted to organize it into some database, but Which one will be effective to store and query this kind of data? Until now I was working only with relational data bases (MySQL), but I guess that is not the solution here (growing number of samples and unique features). Or is it?
Thanks for suggestions!

Not sure what your goal or application is, but if a large array is not feasible due to memory constraints because you're storing mostly zeroes, there's always sparse matrices (https://docs.scipy.org/doc/scipy/reference/sparse.html)

Thanks for the comment. The thing is that I need it as client/server database (that is why I was initially thinking of MySQL). To store bunch of metadata for each sample (and for that MySQL is OK), but also relative abundance. In that case large array is not so much problem of a memory, but table organization, because table structure will have to be very dynamic (constantly increasing number of features, if features are variables).


text:
How do I install hammock?
I'm working with phage peptides libraries. Currently I'm analysing a massive data set from phage selection, NGS readout, preprocessed already (fasta file). I want to analyse the sequences to cluster peptides according to their consensus motifs. I found that Hammock may be helpful to do this.
Here is the GitHub repo for Hammock: https://github.com/krejciadam/hammock. I downloaded the master, unpacked it, but I don't know how to use it.
I guess I need to compile, but I have no idea what files need to be compiled, how to do it, what software I need use to compile etc.
I work with Mac/terminal.

Download a release (Hammock_v_1.2.0.7z) from https://github.com/krejciadam/hammock/releases
Unzip
To run the first example:

cd /your_path/Hammock_v_1.2.0/examples/MUSI/
java -jar ../../dist/Hammock.jar full -i musi.fa

As mentioned in the comments, it's in the manual and any further issues need more information (errors etc).

Answer from @ram-rs, converted from comment:
There's a Hammock.jar file - there' no need to compile Hammock, although they recommend compiling dependencies. Have you tried running java -jar dist/Hammock.jar --help? Follow their steps from Example 1 - it should work.


text:
Determining what RNAseq data is filtered on volcano plot
I am using RNA seq data to analyze genes via a volcano plot comparing differential gene expression of bacteria with and without antibiotic in R. After having created my plot, I am unsure why some of my values which read "0" TPM in one condition and a TPM value that's not "0" in the other condition were not determined to be differentially expressed. Some of the genes on my volcano plot have this difference in TPM and show up as significant on my plot while others with this "0" value difference are not considered to be significant according to my plot.
Sample of my data:
(UO1_D712_##### represents the locus number, the top column represents the different replicates where "1" is untreated and "2" is treated with antibiotic, the letters next to the numbers represent a different technical replicate):
structure(list(`1` = c("U01_D712_00001",

"U01_D712_00002", "U01_D712_00003",
"U01_D712_00004", "U01_D712_00005", "U01_D712_00006", "U01_D712_00007",
"U01_D712_00008", "U01_D712_00009", "U01_D712_00010", "U01_D712_00011",
"U01_D712_00012", "U01_D712_00013", "U01_D712_00014", "U01_D712_00015",
"U01_D712_00016", "U01_D712_00017", "U01_D712_00018", "U01_D712_00019",
"U01_D712_00020", "U01_D712_00021", "U01_D712_00022", "U01_D712_00023",
"U01_D712_00024", "U01_D712_00025", "U01_D712_00026", "U01_D712_00027",
"U01_D712_00028", "U01_D712_00029", "U01_D712_00030"), `1a` = c("6.590456502",
"6.741694688", "7.110342585", "6.299482433", "2.77173982", "2.470330508",
"3.827125186", "6.267229842", "5.524708663", "2.657913228", "2.87209859",
"7.479820548", "4.980185572", "4.210955199", "0", "4.492492822",
"3.611091371", "3.714270433", "7.455036914", "7.045203025", "6.061860857",
"2.925268313", "6.544077039", "5.747318013", "9.97083444", "9.22523089",
"4.20205383", "6.097040679", "2.621192351", "0"), `1b` = c("6.544454427",
"6.601489488", "7.134224619", "5.814043553", "3.280958379", "2.649180803",
"3.860542083", "6.256648363", "5.380427766", "2.581027705", "3.016132165",
"7.405329447", "4.701503289", "4.073814818", "0", "4.304196924",
"3.515977329", "3.843535649", "7.342972625", "6.966606769", "6.122878624",
"3.007522306", "6.495797641", "5.965431621", "9.828050269", "9.219563915",
"4.065778989", "6.105331066", "3.061209408", "0"), `1c` = c("6.608196006",
"6.743010138", "7.102600793", "6.146518601", "3.555184202", "2.364971542",
"4.034053983", "6.158523627", "5.656051812", "2.658660735", "3.054717455",
"7.392164473", "4.950953264", "4.277770477", "0", "4.507936666",
"3.794842979", "3.610794578", "7.471646548", "7.104792624", "6.484767016",
"3.071205184", "6.584425715", "6.20466015", "9.986342122", "9.282943758",
"4.179958213", "6.219551653", "2.984738345", "0"), `1d` = c("6.547155382",
"6.558328892", "6.992501615", "6.449558793", "3.059801464", "2.418800257",
"3.96498952", "6.013208538", "5.279919645", "2.893295085", "1.750510471",
"7.408735671", "4.9425624", "3.804549986", "0", "4.174835979",
"3.806006888", "3.570390524", "7.641137006", "6.976672494", "6.363030106",
"3.083061726", "6.300910093", "6.007490342", "9.926316442", "9.09671588",
"4.320556917", "6.153860107", "2.877230446", "0"), `1e` = c("6.626724417",
"6.577345176", "7.156821278", "6.296873411", "3.618089702", "2.394444986",
"4.129376392", "6.011715246", "5.31197869", "2.00754706", "2.695493528",
"7.538910448", "4.606060035", "3.909472643", "0", "4.346616047",
"3.468681284", "3.338231445", "7.559599613", "7.1527452", "6.232923513",
"3.108624209", "6.535435309", "6.12922864", "10.13108497", "9.310331313",
"3.959568571", "6.182335537", "2.902736258", "0"), `1f` = c("6.419219179",
"6.650459302", "7.319125725", "5.570869357", "2.962845933", "2.55903176",
"4.087597573", "5.995610538", "5.386268651", "2.750800859", "2.416572678",
"7.579148955", "3.952633067", "3.615227674", "0", "4.562838935",
"3.76942104", "3.920096905", "7.935320749", "7.501470652", "6.13700099",
"3.123910608", "6.035971952", "5.706235015", "9.254751395", "8.379630979",
"4.51391973", "5.6890651", "2.43285316", "0"), `1g` = c("6.553221221",
"6.633949852", "7.182305386", "5.769365973", "2.721354972", "1.668390466",
"4.148367057", "6.240883382", "5.458877133", "1.733842637", "2.723803355",
"7.522249899", "4.149567197", "3.780763096", "0", "4.496306813",
"3.645643535", "3.851001768", "7.678552875", "7.283411279", "6.591585956",
"2.879345378", "6.389427003", "5.911222165", "9.851084493", "9.084575304",
"4.272587776", "5.974762147", "2.98852705", "0"), `2a` = c("9.769887737",
"4.550226652", "1.869021464", "4.944848987", "7.9678549", "2.682865013",
"8.495575559", "2.234521659", "3.667196316", "9.180445037", "3.210107621",
"7.21523691", "5.714579923", "5.423986751", "9.118981459", "8.635701597",
"2.742889473", "3.712618983", "8.006057144", "4.999541279", "10.54351774",
"5.880978085", "7.145433526", "6.982416661", "9.339651188", "5.360835327",
"4.699680905", "3.423826225", "6.408271885", "5.038170992"),
    `2b` = c("10.26397519", "1.945664005", "2.086763158", "2.800904763",
    "8.583418657", "1.536094563", "9.32057547", "3.10685839",
    "4.224502319", "10.1252842", "1.811175407", "5.714439316",
    "6.039142559", "3.833361174", "9.757360286", "9.565906731",
    "1.523640473", "2.315033488", "6.312524363", "4.889986456",
    "10.23020108", "4.848727685", "7.533256999", "7.138160378",
    "10.30380331", "3.955469283", "3.167940742", "3.599655687",
    "4.828945262", "3.701043054"), `2c` = c("9.478481216", "2.789289131",
    "3.393949527", "3.754810933", "8.710777154", "1.806170784",
    "9.150005253", "2.612275457", "4.961073313", "9.802701699",
    "2.933183115", "6.532384958", "6.919449225", "4.432699799",
    "9.715063475", "9.265691356", "1.412064593", "3.330131873",
    "6.665979896", "5.158526421", "9.417365584", "4.899531204",
    "8.173459354", "7.271400938", "9.813068613", "4.384622077",
    "3.700645365", "4.457874829", "5.649440022", "3.531010379"
    ), `2d` = c("8.199795497", "2.565711524", "2.202287889",
    "3.856354444", "7.380224849", "2.192476466", "8.14446837",
    "1.144258862", "3.31447122", "8.713146629", "2.697890381",
    "6.304428859", "5.745291803", "4.898396114", "9.173362747",
    "8.339933849", "3.159678152", "4.094587234", "7.608649692",
    "5.280206424", "10.34630403", "5.098585806", "7.262400625",
    "6.150190905", "9.316698845", "5.073027993", "4.695003229",
    "2.485847024", "5.545300465", "4.350571411"), `2e` = c("8.74935033",
    "3.739489484", "0.217205045", "4.413657999", "8.745525588",
    "3.657060086", "9.279834921", "2.898951179", "4.282874018",
    "9.610485827", "3.561102455", "7.228334332", "6.388491443",
    "4.389908652", "8.781086564", "9.178866581", "2.374603596",
    "3.961037408", "7.864369809", "3.654728044", "10.15284858",
    "5.894439123", "7.68020282", "7.12243523", "9.998637438",
    "5.092174395", "4.111530392", "3.776835632", "5.624523213",
    "4.095011377"), `2f` = c("7.648310926", "4.345557215", "1.986576876",
    "4.99426288", "7.087937177", "2.810917253", "7.77706637",
    "2.62822773", "3.581811188", "8.470225989", "3.335757437",
    "7.416094847", "5.208841128", "5.536128034", "8.255571138",
    "7.993319997", "1.9209089", "3.573861828", "7.318814519",
    "5.233804806", "11.05855833", "6.247720809", "6.673407583",
    "6.029960625", "8.806867591", "5.459208493", "4.001428729",
    "3.609936979", "5.876522973", "4.652839671"), `2g` = c("7.555235468",
    "3.892899549", "1.726443458", "5.304546796", "7.039588042",
    "3.027235295", "7.703852207", "1.753183519", "2.909288568",
    "8.385169315", "3.902707541", "7.523315081", "4.978364017",
    "5.49103181", "8.096218606", "7.944822989", "2.352609608",
    "4.155433517", "7.227355741", "5.532668321", "11.24946953",
    "6.159185473", "6.443203375", "5.931761874", "8.7421732",
    "5.502000205", "4.652883503", "3.458323017", "6.566487449",
    "4.89790353")), row.names = c(NA, 30L), class = "data.frame")

Design matrix:
structure(list(SampleID = c("1a ", "1b", "1c",

"1d", "1e", "1f",
"1g", "2a", "2b", "2c", "2d", "2e", "2f", "2g"), Group = c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), Replicate = c("a",
"b", "c", "d", "e", "f", "g", "a", "b", "c", "d", "e", "f", "g"
)), row.names = c(NA, 14L), class = "data.frame")

Here is how I wrote my code to sample the RNA seq data (note ~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv refers to the the sample of data above and ~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv refers to the above design matrix):
R10LB_0vs0.5 <- read.csv2("~/Documents/VOLCANO/R10LB_0ugvs0.5ug.csv", sep=",", check.names = F)
R10LB_0vs0.5 <- janitor::remove_empty(R10LB_0vs0.5, which = "cols")
R10LB_0vs0.5_design <- read.csv2("~/Documents/VOLCANO/R10LB_0vs0.5_designmatrix.csv", sep=",")

rownames(R10LB_0vs0.5) <- R10LB_0vs0.5$`1`
R10LB_0vs0.5 <- R10LB_0vs0.5[,-1]

#colnames(R10LB_0ugvs0.5ug) <- R10LB_0ugvs0.5ug[1,]
#R10LB_0ugvs0.5ug <- R10LB_0ugvs0.5ug[-1,]

R10LB_0vs0.5 <- data.matrix(R10LB_0vs0.5)

R10LB_0vs0.5 <- DGEList(counts = R10LB_0vs0.5, group = R10LB_0vs0.5_design$Group)
str(R10LB_0vs0.5)

Group <- as.vector(as.character(R10LB_0vs0.5_design$Group))
Replicate <- as.vector((R10LB_0vs0.5_design$Replicate))
R10LB_0vs0.5_designmatrix <- model.matrix(~0+Group+Replicate)

keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

R10LB_0vs0.5_Disp <- estimateDisp(R10LB_0vs0.5_filter, R10LB_0vs0.5_designmatrix)

CONTRASTS <- makeContrasts( Group1vs2 = Group1-Group2,
                            levels = R10LB_0vs0.5_designmatrix)

R10LB_0vs0.5fit <- glmFit(R10LB_0vs0.5_Disp, contrast = CONTRASTS[,1])
R10LB_0vs0.5lrt <- glmLRT(R10LB_0vs0.5fit)
R10LB_0vs0.5TT <- topTags(R10LB_0vs0.5lrt, n=nrow(R10LB_0vs0.5_Disp))
write.csv(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.csv")
saveRDS(R10LB_0vs0.5TT, file = "R10LB_0vs0.5_comparison.RDS")

Volcano plot:
ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = PValue)) + geom_point()

p4 <- ggplot(data = R10LB_0vs0.5_comparison, aes(x = logFC, y = -log10(PValue))) + geom_point() + theme_minimal()

#Add vertical lines for logFC thresholds, and one horizontal line for the p-value threshold
p5 <- p4 + geom_vline(xintercept = c(-0.6, 0.6), col = "red") + geom_hline(yintercept = -log10(0.05), col = "red")

# add a column of NAs
R10LB_0vs0.5_comparison$diffexpressed <- "NO"
# if logFC > 0.6 and PValue < 0.05, set as "UP"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC > 0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "UP"
# if logFC < -0.6 and PValue < 0.05, set as "DOWN"
R10LB_0vs0.5_comparison$diffexpressed[R10LB_0vs0.5_comparison$logFC < -0.6 & R10LB_0vs0.5_comparison$PValue < 0.05] <- "DOWN"

p4 <- ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed)) + geom_point() + theme_minimal()

p5 <- p4 + geom_vline(xintercept=c(-0.6, 0.6), col="red") + geom_hline(yintercept=-log10(0.05), col="red")

p6 <- p5 + scale_color_manual(values=c("blue", "black", "red"))
p6 <- p5 + scale_colour_manual(values = mycolors)

R10LB_0vs0.5_comparison$delabel <- NA
R10LB_0vs0.5_comparison$delabel[R10LB_0vs0.5_comparison$diffexpressed != "NO"] <- R10LB_0vs0.5_comparison$gene_symbol[R10LB_0vs0.5_comparison$diffexpressed != "NO"]

ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) + geom_point() + theme_minimal() +geom_text()

library(ggrepel)
ggplot(data=R10LB_0vs0.5_comparison, aes(x=logFC, y=-log10(PValue), col=diffexpressed, label=delabel)) +geom_point() + theme_minimal() +geom_text_repel() + scale_color_manual(values=c("blue", "black", "red")) + geom_vline(xintercept=c(-0.6, 0.6), col="red") +geom_hline(yintercept=-log10(0.05), col="red")

This is the volcano plot I currently have with all of my data:

Can you explain how some of these functions filtered my data (i.e. filterByExpr, estimateDisp, makeContrasts, glmFit, glmLRT)? Why are some of my data points that change from 0 TPM in one condition to some value in another condition showing up on my plot while others are not? Or should I accept that my data analysis is valid or "correct"?

I heard that low expression genes are often filtered out in a first past due to a high chance of false positives, so they could be dropped from filterByExpr. High variance within a condition might undermine the statistical significance of the mean difference across conditions. Is there a way to look at the variance of my  data, and is there a threshold for variance between my replicates that causes them to be filtered out by filterByExpr or other functions?

It is a long time that I don't use R. However, I had a sinilar problem once. Using not cartesian coordinates, ggplot cut several points out of the plot. If you plan to continue using R, search about this problem. Otherwise, if you can accede to it go with MATLAB that has the best tool for volcano plot, by my opinion. It never fails.

edgeR functions do not do any filtering, i.e., they do not remove any genes from the results or from the plots. The only filtering was that which you did explicitly yourself by:
keep <- filterByExpr(R10LB_0vs0.5)
R10LB_0vs0.5_filter <- R10LB_0vs0.5[keep, keep.lib.sizes = FALSE]

For your data, filterByExpr will keep any gene that is expressed at a minimal level in at least 7 samples.
You could relax that requirement, to allow genes to be expressed in only 5 or 6 samples say, by reducing the large.n argument to filterByExpr.

[disclaimer: I understand DESeq2, not EdgeR]
Most likely, low expression values (as in raw counts, prior to any normalisation such as TPM). Expression that has high dispersion across all samples is less likely to be biologically significant.
It's easier to work out what's going on with an MA plot, with average gene expression on the X axis, and log2 fold change on the Y axis.
I'm used to DESeq2 results, so I'd guess the construction would be something like this (assuming that the 'baseMean' field exists in EdgeR results):
R10LB_0vs0.5lrt %>%
    ggplot() +
    aes(y = logFC, col = padj, x=log(baseMean)) +
    geom_point()

It's possible that EdgeR doesn't have anything like baseMean in its default results. As a proxy for this, you could use the count totals from the input matrix (e.g. R10LB_0vs0.5 %>% rowSums() %>% log2()).

Do you mean that the TPM values are too low? Because I believe the ones that have been filtered out have a similar value to other "0" tpm genes that were significant. Also if I use the average of two conditions, will the genes I pull out be less significant because the MA plot will not account for variance between replicates? Thanks for the answer


text:
Question about public availability of human SNP dataset with country of origin
Is there a publicly available data set for humans SNPs, preferably together with the country of origin? If so, could someone please point me towards it.
I've looked through other questions on the site and have been unable to find a similar one. I don't know how sensitive human SNP data is, so I'm not sure if this would be public.

1000 genomes, HGDP and SGDP all contain SNP data from individuals where the country of origin is known.
They all have slightly different numbers of people from different populations (1000 genomes and HGDP have more individuals in each pop, but fewer pops, and SGDP is the opposite).

Thank you, this is very helpful. Do you know if any of these websites have .vcf files already ready to download? Or will I need to follow a pipeline to convert the .fastq files to .vcf?

Yes the 1000 genomes files can be found on this ftp mirror http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/data_collections/1000G_2504_high_coverage/working/20201028_3202_raw_GT_with_annot/ - although the total vcfs tend to be quite big

hgdp vcf can be found here ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/


text:
Converting FASTA Sequences to SNP FASTA Format
I have two FASTA sequences. Each is from a different species' mitochondrial genome. The sequences are pasted below.
>speciesA_control_region_partial_sequence_mitochondrial
TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA
>speciesB_control_region_partial_sequence_mitochondrial
TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA

My goal is to extract all SNPs in these sequences and format them in SNP FASTA format, as shown below:
>rs11540652
GAACC[A/G]GAGGCCC

I specifically need this format (rather than, say, VCF) so I can search for RFLP sites using the following website:
http://bio.kuas.edu.tw/snp-rflping2/rflpUI.jsp
Does anyone know of an existing program which will accomplish this? Also, what additional materials do I need to detect these SNPs? I presume that I need assembled mitochondrial reference genomes for each species - anything else?

Do you have a VCF file? It should be possible to compute one using bcftools and it would give the position in the reference sequence so generating these “SNP” sequences could be done using this and the reference sequence. The full sequence for both samples will take longer to compute and is not necessary.

Your question is not very clear. Where does the flanking sequence in :
>rs11540652
GAACC[A/G]GAGGCCC

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the [] is from the reference... I don't know an existing existing program but here's a bespoke one.
speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

Gives:
>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.
I hope that helps.


text:
How is the instantaneous rate of an instantaneous rate matrix determined?
In molecular evolution, I know that the instantaneous rate of the rate matrix is the limit of the rate as time approaches 0, but what I don't understand is how this rate is established or what I can do with this number.
For instance, if a discrete character is thought to have changed from 0 to 1 five times on a tree with a root age of 40 million years, is the proportion 5 / 40 million the same as the instantaneous rate of change? If not, how exactly are they different?
From a different perspective, let's assume an (absolute) instantaneous rate of 0.0001 and let's say I have a tree with a root age of 2 million years. Can I multiply 0.0001 by two million to get the expected number of transitions of that character across the tree?

related: https://bioinformatics.stackexchange.com/questions/13605/what-does-instantaneous-rate-matrix-mean

@MaximilianPress I saw this post. Although both questions concern rate matrices, they are nevertheless quite distinct.

I agree; note that I didn't mark it as duplicate. However, they have a substantial overlap and it is worth annotating as such, for anyone coming across this question.


text:
plink doesn't generate a sexcheck-file
I was executing the following command plink.exe --file GWAS --check-sex --noweb
Expecthing to get a plink.sexcheck file but I don't get anything. Not even an error. The output is always:
Skipping web check... [ --noweb ]
Writing this text to log file [ plink.log ]
Analysis started: Sat Mar 13 12:24:18 2021

Options in effect:
        --file GWAS
        --check-sex
        --noweb

75 (of 75) markers to be included from [ GWAS.map ]

Do you know what could be the issue there? I have a GWAS.map and GWAS.ped file. Do I need any additional information for this to work?


text:
Biopython BioSQL error : AttributeError: 'DBSeq' object has no attribute '_data'
Hi I am playing with Biopython and BioSQL,
I manage to create an MySQL database using Managing local biological databases with the BioSQL module tutorial:
and load 2909 fasta sequences in it through:
db = server["test_1"]

count = db.load(SeqIO.parse("Newtextfile-03.txt", "fasta"))
        print("Loaded %i records" % count)
        server.commit()

from this file: Newtextfile-03.txt
problem is trying just to print the database index list with:
for i in db:
    print(i)

works , while trying print element by element:
for i in db:
    print(i, db[i])

or with:
for key, record in db.items():
            print ("Key %r maps to a sequence record  %s" % (key, record))

gives an error:
2557 Traceback (most recent call last):
  File "./try_001.py", line 142, in <module>
    f()
  File "./try_001.py", line 120, in f
    print(i, db[i])
  File "......../lib/python3.8/site-packages/Bio/SeqRecord.py", line 652, in __str__
    lines.append(repr(self.seq))
  File "........//lib/python3.8/site-packages/Bio/Seq.py", line 109, in __repr__
    return f"{self.__class__.__name__}({self._data!r})"

AttributeError: 'DBSeq' object has no attribute '_data'

the error occurs for three records namely : [2164, 2522, 2557]
inspecting the Newtextfile-03.txt file results in finding fasta sequences of less than
100 aminoacids.
Is there any way to resolve the error ? It sounds strange to me that I can load the records into the database but not retrieving them. Is anything I am missing about the record load ? or similar ?
I am working with the database using:
 server = BioSeqDatabase.open_database(
    driver='MySQLdb' ,
    user=password.userz,
    password=password.pazzword,
    host="XXXXXXXX",
    db="bio-data-first")

db = server.new_database("test_1", description="Just for testing")
        server.commit()

using  driver="MySQLdb", wasn't able to have driver='mysql.connector' to work,
biopython is version 1.78

This does look like a bug but the code has been refactored a lot recently. Try running with biopython 1.79 (it's not released yet so you'll need to install directly from the github)

so its about 'DBSeq' object ?

thanks @Chris_Rands I'll give it a try as soon as I figure out about virtual enviroments and how to install modules from code

found out that: for key, record in db.items():
print(key)
seq_record = db.lookup(gi=record.id)
print(seq_record.seq) works for all the elements

It's a bug in Biopython 1.78.
What those three sequences have in common is that they're less than 61 characters long.  When you print the record it implicitly calls __str__ for the record, which calls __repr__ for the Seq.  DBSeq gets its method from Seq, which looks like:
    def __repr__(self):
        """Return (truncated) representation of the sequence for debugging."""
        if len(self) > 60:
            # Shows the last three letters as it is often useful to see if
            # there is a stop codon at the end of a sequence.
            # Note total length is 54+3+3=60
            return f"{self.__class__.__name__}('{str(self[:54])}...{str(self[-3:])}')"
        else:
            return f"{self.__class__.__name__}({self._data!r})"

...and then it only tries to access the nonexistent _data for the short sequences.
So there's a mismatch between DBSeq's lack of a _data attribute and its use of Seq.__repr__.  I don't see a DBSeq class at all in the latest Biopython, so maybe this is fixed already.  But there's no tag for 1.79 yet that I can see and I can't get the very latest code to import successfully so I'm not sure.  In the meantime if you avoid letting repr(record) happen you should be able to work around it.


text:
Doubt about using TPM for statistics
We designed an experiment to explore the potential role of carbon dioxide on algae physiology using RNA-Seq. We analyse the differential gene expression using DESeq2 but now we are interested into analyse a few genes in detail.
My PI say to make the statistics over the ration response (expression_treatment/expression_control - 1) and my question is the following: Can I use TPM to calculate this ratio and make following statistics? I read in the documentation of DESeq2 that TPM we can't use this normalisation method to differential expression and my my other question is Why shouldn't?
My advisor ask to re-do the statistical analysis using the ratio response (rr) so i need the rr for each sample. I think this doesnt make sense cause the result will be or should be the same that using DESeq2.
Update: I'm using the fold-change calculated by DESeq2 and I'm not sure if I have transformed the data correctly. I took the two columns log2FoldChange and lfcSE (both are log-transformed variables) so I went back to the non-transformed variable as follow: FC=2^log2FoldChange and SE=lfcSE*2^log2FoldChange. I saw this approximation in this post. Is that correct? I have obtained a very broad range of error/

This is a common question, and the answer is while you can calculate this, it's not statistically robust.
You're likely to arrive at false conclusions.
Why is this?
Because while TPM is something you can calculate from a single sample, your question involves multiple replicates and those replicates vary between each other.
Counts in units of TPM are missing information about the variation between samples.
This extra information affects fold change estimates, and is why differential analysis of sequencing data is complicated.
See this answer to a related question, and this often-cited blog post for more information.


text:
In SBOL, is there a way to distinguish biological replicates?
I have a library of biological implementations (i.e. plasmids) of genetic circuits and their experimental data associated. I want to implement an in silico version with SBOL, but to be consistent with what I have I'm my freezer, I need to identify not only the genetic circuit but also the biological replicate (i.e. the different colonies picked), as they have their own associated experimental data.

Support for distinguishing replicates was added in SBOL 2.2, along with other support for linking data throughout the experiment cycle.  I will explain this in terms of SBOL 3, which is simpler, then give the SBOL 2 translation at the end.
For SBOL 3, all designs, including those for a plasmid, a strain, etc, are specified in terms of the Component class. This is just a "paper" design, however, and may or may not actually exist. When you have a real physical aliquot that you want to talk about, this is represented by an Implementation object: this is your biological replicate.
The Implementation doesn't really have any content of its own; its nature is defined by linking back to the intended design with a prov:wasDerivedFrom link. If you've verified that what you got is what you intended, you can supplement that with a built link as well.
One additional note, however: when you're representing a colony, you're probably not actually representing just a plasmid, but rather are representing a particular strain that has been transformed by a plasmid. The actual data you have may then be one level further removed, since the data has been gathered from a sample under particular growth conditions. SBOL can represent that as well, by describing the sample as a mixture of particular media, inducers, etc., with a transformed strain, which is itself a cell type containing your plasmid.
Here is an example diagram showing just this sort of relationship, of data associated with three replicates of a sample.

You can find this example in the SBOL3 tutorial from COMBINE 2020.
Finally, in SBOL 2, this works exactly the same way, except that instead of a Component your implementation points to a ComponentDefinition or a ModuleDefinition, depending on the nature of the design you're representing. I won't go into that distinction, because it was a bad idea and there's a reason we got rid of it for SBOL 3.

Indeed, as @jakebeal says, the "colony" is a strain transformed with a plasmid, and then this Implementation can be used in experiments with different media and inducers.

Exactly, @AlejandroVignoni: you can end up with trees of replicate derivation. I personally value that quite a bit, because the "flat" notion of "biological replicate" vs "technical replicate" often ends up obscuring these relationships.


text:
Is there a good practice (or an easy one for a SBOL newbie) to set up a collection of biological parts and devices?
I have a collection of biological parts and devices in a particular format that basically stores information of the sequence of the elements. My intention is to move to a more standardized data structure with a particular interest in Synthetic Biology. SBOL 3 looks like a good choice, then the question is if there is a way to kind of "import" data from my sequence-based database into SBOL?

For SBOL 2, you can readily import data via the online converter tool at https://converter.sbolstandard.org/
The converter basically just wraps the conversion function found in the SBOL libraries, so if you're coding, you can do the conversion directly from the libraries.
For SBOL 3, there isn't yet an online tool, but James McLaughlin has a converter in his JavaScript library: https://github.com/sboltools/sbolgraph
I believe he is working to bring it online shortly, so hopefully this answer will soon need to be updated or replaced by an answer from him.
An anecdote on the difficulty level: I recently started using an SBOL3-->SBOL2-->GenBank conversion using pySBOL3 and pySBOL3. It took a couple of hours to setup the pipeline, and isn't perfect, but it's pretty doable (and ought to be generalized for others to use as well at some point). The other way would be easier, since GenBank contents are pretty simple.

Further to jakebeal's answer there is also an excel2sbol library (still in alpha) that can take input to a standard template and convert to SBOL2.
For SBOL3 there is also a python library, but using this to create SBOL3 is more involved.


text:
How do you find differences when your metagenome sample size is small?
I received a matrix of genomics features (KEGG annotated metagenomes) abundance from 6 samples belonging to 2 groups. My collaborator is interested in finding differences between the two groups.
The matrix, after scaling and removing zero-variance, is roughly 4000 features. Via shapiro test I already know some 30% of the features are not normally distributed across the samples, so I assume PCA is out of option and any test with the normal distribution assumption.
I already made some preliminary analysis using NMDS to cluster my samples and was thinking of trying also t-SNE and recursive feature elimination. However I am not sure how to find "significant" differences given the small sample size. Any suggestion? Would something like LEfSe work in this instance?

I think there are several ways to figure this out. I wouldn't worry too much about normality as most of the methods are either nonparametric or pretty robust to violations.

Since you are using KEGG, I would suggest working at the pathway or module level of annotation (if your features are e.g. KEGG Orthology groups). In principle, you can then combine information across features as well as across samples. Simplest I can think of, without even looking at the literature, would be to use permutations (shuffle row and column labels) to look for modules different in abundance between groups. (Many existing methods are doing some more sophisticated version of this.)
I think that LEfSe in theory should work quite well here, it's been a while since I thought about it but I believe this is the problem that it is designed to solve, assuming that I understand your question.
There are a lot of methods that are previously described for doing this kind of work. I would suggest reading some reviews of methods, for example the section "Multivariate differential abundance testing" in this article, or here, looking for some inspiration. To some extent it depends on which questions you are interested in answering and what exactly your features represent.

Calle's article has been greatly helpful! I did not know of all the work done on compositional data analysis. Now I have plenty of new tests to try. Thanks.


text:
How to interpret track height in Chipseq?
I recently received some .bigwig files from a chipseq experiment.
I loaded the files into IGV.
The interpretation is rather intuitive.
My question is the following:
What exactly do the track heights represent?
For IgG control, the track height is low ranging from 0~5.
For positive controls, such as H3K4me3, I see clear peaks with height ~100.
Do track heights simply mean the number of reads mapped to that genomic interval during peak calling?
Does it literally mean there are ~100 reads that aligned to that region in the H3K4me3 .bam file?
Thanks so much.

Depends how the files were made. In the simplest case yes, the height represents the pipeup of reads from the BAM file that was used. For a direct (visual) comparison you have to normalize the files though as otherwise sequencing depth confounds the height of the peaks, like if file A has ten times more reads and the peak is 10 times higher than in B then this is probably only due to differences in total read number (=sequencing depth). I wrote down my preferred method to do that in this biostars post: https://www.biostars.org/p/413626/#414440
Does this make sense to you?


text:
Common gene naming conventions and how to convert between them
I have two questions, one is a direct question about some RNA expression data.  The other is more broad, but motivated by the first.  I have downloaded an expression set for a class project and I am unfamiliar with the naming convention, I recognize entrez ID but not this one.
    > rownames(Rao)[1:10]
 [1] "0610007P14Rik" "0610009B22Rik" "0610009D07Rik" "0610009L18Rik" "0610009O20Rik" "0610010B08Rik" "0610010F05Rik"
 [8] "0610010K14Rik" "0610011F06Rik" "0610012G03Rik"

So my first question is what is this gene naming convention.  The second is what are many other conventions I should recognize and what tools can I use for converting between them?  Lastly what kind of errors may crop up when I do that and how can I check for mismappings or non-unique mappings?
I am primarily interested in R but methods in python would also be appreciated!  Its always good to have more tools.

These look like Mouse Gene identifiiers. You might be able to use org.Mm.eg.db to map them to aliases. There should be a way to use biomaRt as well, but I don't know the speciifics - maybe a mus musculus gene symbol mart is available.

Thanks for reminding me to check biomaRt.  I think it has what I need, just reading a blog post about it now.  https://yiweiniu.github.io/blog/2018/09/Biological-ID-Conversion/

Biologists are familiar with SYMBOL (for good reasons).
Bioinformaticians like ENTREZID ENSEMBL REFSEQ GENBANK etc, I'm definitely missing some here.
You don't need to know these by heart, just know that databases don't map 1:1.
I've worked with mouse RNAseq data before, so I agree with user::RamRS comment that it looks like a mouse gene. What type of identifier is it? I had no idea.
You can find this out at NCBI::GENE, by querying for, let's say Rao[1] "0610007P14Rik". (https://www.ncbi.nlm.nih.gov/gene/?term=0610007P14Rik)
NCBI::GENE tells you the SYMBOL::Erg28, ENTREZID::8520, and in the far right a column called ALIAS::0610007P14Rik.
Now that you know your identifier type (ALIAS), you can map between databases in R.
There are many packages to do this. Someone mentioned biomart in the comment.
I use AnnotationDbi
library(AnnotationDbi)
library(org.Mm.eg.db)

Rao <- c("0610007P14Rik", "0610009B22Rik", "0610009D07Rik", "0610009L18Rik", "0610009O20Rik", "0610010B08Rik", "0610010F05Rik", "0610010K14Rik", "0610011F06Rik", "0610012G03Rik")

#use the function select to retrieve the column type "SYMBOL" with keys "RAO".
#"RAO" is a keytype of "ALIAS". And I want to search in the mice database "org.Mm.eg.db"

AnnotationDbi::select(x = org.Mm.eg.db,
+                       keys = Rao,
+                       keytype = "ALIAS",
+                       columns = "SYMBOL")
'select()' returned 1:1 mapping between keys and columns
           ALIAS        SYMBOL
1  0610007P14Rik         Erg28
2  0610009B22Rik 0610009B22Rik
3  0610009D07Rik         Sf3b6
4  0610009L18Rik 0610009L18Rik
5  0610009O20Rik         Dele1
6  0610010B08Rik        Gm2026
7  0610010F05Rik 0610010F05Rik
8  0610010K14Rik 0610010K14Rik
9  0610011F06Rik       Mettl26
10 0610012G03Rik 0610012G03Rik

In the cases when they don't map 1:1, the function will throw a warning, and you can access it by warnings(). It is rarely important though.
Hope this answers the question.

Thanks that does help quite a bit.  Is there a way to auto remove all the mismatches?

Like with some kind of binary vector so i can filter out the mismatches?

pipe it to dplyr::filter or subset and use the criterion: ALIAS != SYMBOL


text:
What is the best way to programmatically convert Ensembl ids from one release to other? e.g. transcript ids from release 75 to 100
3 solutions that I had to unfortunately rule out.
Potential solution #1
As far as I know, the conventional way is to do this is by using CCDS ids (e.g. using biomart). However, the issue is that, as I see on the CCDS project page, the last release of CCDS ids was made almost 3 years ago (!) and it seems to be compatible with Ensembl release 92. Then how should one convert ids to the latest Ensembl release e.g. 100?

CCDS Release
Species
NCBI Annotation Release
Ensembl Annotation Release
Assembly Name
Assembly ID
Made Public

22
Homo sapiens
109
92
GRCh38.p12
GCF_000001405.38
06/14/2018

Potential solution #2
I am also aware of Ensembl's ID History Converter. However, I am looking for a programmatic approach. Also, I noticed that if I try to convert ~>20K ids with ID History Converter, in my experience, there is always some issue. e.g. an error I got when I tried to convert ~70K ensembl transcript ids.

As a side note, I wonder why this ID conversion functionality is not available via Ensembl's REST API. To me, id conversion seems to be a basic enough task that almost every bioinformatician has to do. So it would be immensely helpful to be someday able to do this programmatically using Ensembl's mighty a REST API.
Potential solution #3
Finally, as another potential solution, over at Ensembl's FTP website, located within the directory of the release to which I am interested in converting my ids to (release 100, target release), I saw that according to a README file, "ensembl_archive_" directory should contain "Data on historical Ensembl releases". So, I wondered if this directory contains a map between ids from different releases. However, without any extra info about column names, contents, file formats etc, I couldn't really figure out how to use those files. So, I had to rule this one out too.
Index of /pub/release-100/mysql/ensembl_archive_100/
../
CHECKSUMS                                          26-Mar-2020 14:31                 132
ens_release.txt.gz                                 26-Mar-2020 14:31                1098
ensembl_archive_100.sql.gz                         26-Mar-2020 14:31                 943
release_species.txt.gz                             26-Mar-2020 14:31               39620
species.txt.gz                                     26-Mar-2020 14:31                6669

Of note, I saw a similar "mysql" directory over at the directory of the release whose ids I want to convert (source release): http://ftp.ensembl.org/pub/grch37/release-100/mysql/homo_sapiens_core_100_37/ . However, again, I do not understand its contents without any extra information.
So, finally, I am not sure what is the best way to programmatically convert Ensembl ids from between different releases. As a test case scenario, I want to convert all the protein_coding transcript ids from GRCh37 (release 75) to GRCh38 (release 100). If anybody has any suggestions about this, please kindly let me know.

The ID history converter is available as a Perl script.
This script accesses these MySQL tables, so you could query these tables directly.
Note that the mapping session is only one release at a time, so you would need to run 25 queries to go from release 75 to 100.

Thank you for your answer. To use the perl script, I wonder if there is any tutorial to refer to. I am particularly interested in knowing how to "tell" the script to convert from a given specific release to other specific release e.g. from 75 to 76.

It just gives you the whole history, you would need to parse out the releases you want.

IDmapper.pl does the job, but I wish there was an easy way to install its dependency called Bio::EnsEMBL::Registry. I feel that perhaps a wrapper of sorts, for example in the Anaconda system, would be immensely helpful. (Please somebody make it happen! :D )

Here's one way to do it programmatically using the biomaRt package (your potential solution #1):
library(biomaRt)
library(dplyr)

mart37 <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                  host    = "grch37.ensembl.org",
                  path    = "/biomart/martservice",
                  dataset = "hsapiens_gene_ensembl")

mart38 <- useMart(biomart = "ENSEMBL_MART_ENSEMBL",
                  path    = "/biomart/martservice",
                  dataset = "hsapiens_gene_ensembl")

#pick all genes on chr22 from build 37 as a test dataset::
out37 <- getBM(attributes = c('hgnc_symbol', 'chromosome_name',
                              'ensembl_gene_id', 'ensembl_transcript_id', 'ccds'),
                  filters = c('chromosome_name'),
                  values = list("22"),
                  mart = mart37) %>% #empty string CCDS ID's
  filter(ccds != "")

#NOTE:: to check which filters can be used to query mart37, run this command:
#View(listFilters(mart37))

out38 <- getBM(attributes = c('hgnc_symbol', 'chromosome_name',
                              'ensembl_gene_id', 'ensembl_transcript_id', 'ccds'),
                        filters = c('ccds'),
                        values = list(out37$ccds),
                        mart = mart38)

What this does is query the build 38 DB via the CCDS ID's of build 37. You can include other variables in this out38 dataframe by simply adding to the attributes argument above. And in order to see all of the potential variables that can be included, you just have to run View(listFilters(mart37)) or View(listFilters(mart38)), depending on which build you are interested in.
Does this solve your problem?

your solution could indeed be useful. I am just curious which ensembl release numbers are used when mart37 and mart38 were built. Is there a way to specify that?

mart37 corresponds to the GRCh37 release, while mart38 corresponds to GRCh38. This is specified in the host= argument in the code above. mart38 is the default, and so is not specified. Does this answer your question?

P.S. if you're referring to the version number of the ensembl release, you can obtain this by adding 'ensembl_gene_id_version' to the attributes array in the code above.

By Ensembl release number I meant version numbers of the assemblies maintained by ensembl, e.g. release 75 corresponds to GRCh37 for human genome.

ensembl_gene_id_version could be useful if it is possible to know to which release that belongs to. e.g. how would I know if x.3 id belongs to genome release 90 or 91? To get this information, I would have to run 'Ensembl IDmapper' which I would like to avoid, because of the issue I mentioned in my question.


text:
How to check what software has been used to generate bam file?
I have several bam files from unknown origin. I'd like to know what software was used to obtain those bam files, e.g. bwa-mem. Is it possible to somehow check it?

If you're lucky the information is stored in the BAM header in the lines starting with @PG which stores the command lines run on the BAM, use samtools view -H your.bam to display the header.

Indeed, there was the full command line, thanks!


text:
Nextflow error: can not invoke method minus () on null object
Hi I was trying to use plink in nextflow and getting the error mentioned in the subject. Here is my code, any help what I am missing. I am new in nextflow, it was ok when input was a single file, but for grouped files, I am getting lost.
params.dir='data/'
params.publish="output/"
Channel
.fromFilePairs("${params.dir}*.{bed,fam,bim}",size:3){file -> file.baseName}
.set {plink-data}

process test1{
publishDir params.publish
input:
set t1,file(plk1) from plink-data
output:
file "${t1}.{bed,fam,bim}" into test1_results
script:
output_file = "${t1}"
base = plk1[0].baseName
"""
module load plink
plink --bfile $t1 --maf 0.05 --make-bed --out ${output_file}
"""
}
process test2{
publishDir params.publish
input:
set file(bed), file(fam), file(bim) from test1_results.collect()
output:
set t2,file '${plk2}.{bed,fam,bim}' into test2_results
script:
output_file = "{t2}"
base = plk2[0].baseName
"""
plink --bfile $t2 --maf 0.01 --make-bed --out ${output_file}
"""
}

Channel names must consist only alphanumeric characters including underscores. Your channel called 'plink-data' contains a hyphen. This is problematic because it is interpreted as a minus() operation to remove 'data' from 'plink'. You get a NullPointerException since 'plink' hasn't been defined. You can test this easily with the following example:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr-names }

chr-names.view()

Results:
Cannot invoke method minus() on null object

The solution is to just use an underscore instead of a hyphen:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr_names }

chr_names.view()

Follow-up issues:
I think if you want to use the bfiles plink option (and I think this is reasonable), I think it makes the most sense to just keep your bfiles in a tuple:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

If you need to reference each bfile separately, you could just unpack the tuple and reference each file accordingly:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

Note that in the first example, bfiles could be replaced with "*" since no reference is made to the tuple of bfiles anyway. Also note that the list of files produced using the fromFilesPairs operator are returned in lexicographic order.

Thanks @Steve, I modified the code but getting a new issue:
'process test1{
publishDir params.publish
input:
set t1, file(plk1) from data
output: set t2, file "${t2}.{bed,fam.bim}" into test1_results
script:
"""
plink --bfile $t1 --maf 0.05 --make-bed --out ${t2}
"""
'
the error is: no such variable: t2
I have 3 files (.bed,.bim,.fam), process test1 will generate another 3 (.bed,.bim,.fam) files which will be used in the next process.
How could I declare t1,t2,t3 etc (the inputs and outputs) correctly?

No worries @zillurrahman, I have added some example code to show how you could declare your inputs and outputs. The inputs could just be a tuple of bfiles or each file could be referenced separately. I haven't used plink before and have only had a (very) quick read of the docs, so there may be some more issues. Hopefully it's enough to get you started.

Thanks a lot @Steve. It works as I expected! One thing, plink also generate a log file, I need this as one of the outputs but the log will not be used in the next process. Now I am adding more processes in the script, I will get back to you if I got stuck again.

So for Nextflow to be able to publish a file using the publishDir directive, it (the file) needs to be referenced in your process' output declaration block somehow. If it's not required by any downstream process, you can just add another output channel, for example: path("${sample_name.log}") into test1_logs. Please let me know if that doesn't make sense. Cheers.

Thank you @Steve, It seems I can declare multiple output, which is very much useful for my pipeline. Just one thing, if I  'path("${sample_name.log}") into test1_logs' it shows 'No such variable: log'.

Erp there's a typo in there.. you'll need: path("${sample_name}.log") into test1_logs if your log's filename fits that format. Whether it does or not, you'll need to check. HTH.

Let us continue this discussion in chat.

Thanks @Steve, I also tried it. If I use path("${sample_name}.log") , it shows error execution prosess; caused by missing output files. If I use path("${sample_name}".log) , it shows no such variable: log

What's the name of your log file? I haven't run plink before, so I just used "${sample_name}.log" as an example. If the name of the log file is 'plink.log' then you'd use: path("plink.log") into test1_logs in your output declaration.

Thanks @Steve, The final output.log is the log file. Sorry I failed to clarify. However, path("${sample_name}-t1.log") into test1_logs works as I expected! Thank a lot again.


text:
Calculate genome coverage and depth from alignment
I have a .bam alignment file and a genome reference .fasta file. I am looking for a easy to use tool (that I can reference in a publication) to calculate the percentage coverage of the reference by the alignment and if possible also: average depth of reads that cover the gene, min/max depths. Running plot would be nice but not essential.
I wrote some python code to calculate coverage percent, but I want to reference the tool I use in a paper and dont want to mess with making a git repo and linking to it in the paper just for this simple analysis.

You could use samtools coverage as explained in the manual of samtoools.
Here is a example which is also described on the manual site.
samtools coverage -r chr1:1M-12M input.bam

#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

Are you sure that those statistics are the ones that make the most sense for RNASeq?

I've used https://genome.sph.umich.edu/wiki/BamUtil before. I think this has all the functionality you want and the Devs are responsive which is nice.


text:
Nextflow: No signature of Method: chr() is applicable for argument types; unknown variable
I was trying to use plink in nextflow. My first process will split the input (single bed|bim|fam file) into 22 chromosomes (each with bed|bim|fam files) which would be used in subsequent processes. My code and error log is below. Any help?
params.dir='data/'
params.publish="output2/"
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}

Channel
.of(1..22)
.set {chr_names}

process step1a {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from data
output:
tuple val("${prefix}"), path("${prefix}.step1a.{bed,bim,fam}") into step1a_results
tuple val("${prefix}"), path("${prefix}.step1b.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.step1b.bim") into step1a_bims
path("${prefix}.log") into step1a_logs
path("${prefix}.step1a.dupvar") into step1a_duplicates
path("${prefix}.step1a.log") into step1a_log
path("${prefix}.step1b.log") into step1_log
script:
prefix=chr"${chr_names}"
"""
plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
"""
}
process step1b {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from step1a_out
tuple val(bims), path("${prefix}.step1b.bim") from step1a_bims
output:
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
path("${sample_name}.step1c.log") into step1b_logs
path("${sample_name}.step1b.snplist.txt") into step1b_snplist
script:
"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
plink --bfile "${prefix}.step1b --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""
}

The error message is:
N E X T F L O W  ~  version 20.10.0
Launching `t4.hf` [zen_picasso] - revision: b0ff8db413
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
  No signature of method: chr() is applicable for argument types: (org.codehaus.groovy.runtime.GStringImpl) values: [DataflowQueue(queue=[DataflowVariable(value=1), DataflowVariable(value=2), DataflowVariable(value=3), DataflowVariable(value=4), DataflowVariable(value=5), DataflowVariable(value=6), DataflowVariable(value=7), DataflowVariable(value=8), DataflowVariable(value=9), DataflowVariable(value=10), DataflowVariable(value=11), DataflowVariable(value=12), DataflowVariable(value=13), DataflowVariable(value=14), DataflowVariable(value=15), DataflowVariable(value=16), DataflowVariable(value=17), DataflowVariable(value=18), DataflowVariable(value=19), DataflowVariable(value=20), DataflowVariable(value=21), DataflowVariable(value=22), DataflowVariable(value=groovyx.gpars.dataflow.operator.PoisonPill@28c22659)])]
Possible solutions: run(), run(), any(), is(java.lang.Object), use([Ljava.lang.Object;), tap(groovy.lang.Closure)

Source block:
  prefix=chr"${chr_names}"
  """
  plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
  plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
  plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
  plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
  """

Tip: when you have fixed the problem you can continue the execution adding the option `-resume` to the run command line

If I use "prefix=chr"${chr_names}" " inside the code block (""") it shows: unknown variable 'prefix'

The documation re conditional scripts is useful here:

Process scripts can contain conditional statements by simply prefixing
the script block with the keyword script:. By doing that the
interpreter will evaluate all the following statements as a code block
that must return the script string to be executed.

But the statements don't necessarily need to be conditional. I.e. they can just be declarations like in your example. You get the error 'No signature of Method: chr()' because Nextflow interprets 'chr' as a method call when you use:
prefix = chr"${chr_name}"

If your intention is to simply add the 'chr' prefix to each chromosome name, you'd need:
prefix = "chr${chr_name}"

Notice that I changed your variable name. This is because in your code, you still need to feed in the channel of chromosome names. The easiest way to do this is to use the each qualifier:
process myprocess {

    input:
    tuple val(sample_name), path(bfiles) from plink_data
    each chr_name from chr_names

    script:
    prefix = "chr${chr_name}"

    """
    echo "My sample name is: ${sample_name}"
    echo "My prefixed chromosome name is: ${prefix}"
    """
}

Unless you need chromosome names without the 'chr' prefix, another solution might be to simply add the prefix to each one initially and simply avoid having to prefix them later:
Channel
    .of(1..22)
    .map { "chr${it}" }
    .set { chr_names }

It seems process 1 and 2 running in random order. Is there any way to specify that after finishing the first process, the second process will start. Because the outputs of the first process will be the input of the second process. In my case my second process shows error because it does not find an input, which is an output of the first process. Did I clarify the matter?

I suspect the problem arises because you've specified two input and output channels ('step1a_out' and 'step1a_bims'). It's not clear to me why you need the two channels, since the 'step1a_bims' looks like a subset of 'step1a_out'? Please see understanding how multiple input channels work

Ohh! Thank you. Just one thing, How could I declare prefix of the previous process output as an input of the current process. I have chr(1-22).step1b.bed|bim|fam as output of process step1a; Now I want to declare only chr(1-22) as input in the step1b process. I want different suffix inside process step1b. I am not sure if I could clarify the issue.

Does my answer to your most recent question answer your question?

Yes. Thank you. Still wondering how could I declare prefix (characters before first .) as input in the next process which is the output of the previous process.

No worries at all @zillurrahman. If I understand correctly, the easy way would be to include it in your output tuple as an extra value for your prefix. But if it's already part of (i.e. the prefix) of another value then you could again use the map operator to get it back. Perhaps something like: map { some_name, some_thing -> tuple( some_name.substring(0, some_name.indexOf('.')), some_name, some_thing ) } ? A concrete example would be great if the above doesn't help.

Thanks a lot @Steve. I tried different ways but failed, maybe because I am new in nextflow. My files are chr1-22.step1b.bed|bim|fam in out_ch1. I need a way to find chr1 to chr22 as prefix for the next input so I can modify the suffixes. for example input chr1-22.step1b.bed|bim|fam --some command --out chr1-22.step1c.bed|bim|fam into out_ch3.
Also how can I use subset of the outputs channel. From above example I need only chr1-22.step1b.bim as a $bim --some command to run -out chr1-22.step1b.snp.txt. So far I tried "filter" but failed. Thanks again for helping all the way.

@Steve, I tried to filter as: tuple val(bim), path(bims) from step1a_bim.filter{bim -> bims=~ /*bim/} But I am getting error unexpected input: Also I am still not getting the prefix correctly. Perhaps any example would be helpful. Do you know any training module to run plink, gwas in nextflow?

Sorry @zillurrahman - I'm not sure I can get you a useful answer without a minimal example of exactly what you're trying to do. On file-like objects you can call getSimpleName() to get the prefix, see check file attributes. The problem is that this won't work on a tuple of files. I guess I'm confused as to why you can't just include the chromosome name in the output tuple of the first process? This way you'd avoid having to try and find the prefix in the first place.

Using my answer above as an example, your output declaration could look like: tuple val(sample_name), val(chr_name), path("${sample_name}.${chr_name}.step1b.{bed,bim,fam}") into step1a_out. Also, I'm not aware of any good plink/gwas pipelines but you could try having a look here: https://github.com/nextflow-io/awesome-nextflow. HTH.

Thanks a lot. In that case how can I declare inputs in the next process? If I use step1a_out two time (first for bfiles, and then for bim) I got error Channel 'step1a_out' has been used twice as an input by process 'step1b'

If I understand correctly, the input declaration for the next process will be basically the same as the first processes' output declaration. If that's the case, you'd need: tuple val(sample_name), val(chr_name), path(bfiles). Since the bfiles tuple will contain the bim  file as the second element, you should be able to access it directly with ${bfiles[1]}. Alternatively, in your script block, you could unpack the bfiles: (bed, bim, fam) = bfiles, then just use ${bim} as your variable. If these variables don't need to be global, consider instead: def (bed, bim, fam) = bfiles. HTH.

Thanks a lot @Steve . It works. Just one thing. How could I specify output directory of my results? I tried params.publish='output2' then in each process publishDir params.publish. I got only link of the output files in output2. How could I get actual result files in my desired directory? Again thanks a lot for helping me all the way.

Ohh! I got it. I just need to mention mode: 'copy' in publishDir declaration. Thanks!


text:
Nextflow: ordering the execution of processes; outputs of process A are inputs of process B
I need to build a pipeline where 2nd process starts after completion of the 1st process and the 3rd after 2nd. As I understand nextflow's processes run parallel. But I need it to run in a order 1st>2nd>3rd processes and so on. Because each process runs on an array and produce multiple outputs which are used by the next. If they run paralleled, the 2nd process may look for input 20 whereas the 1st process produced output 10 and then the process stops, because the command in the 2nd process does not find the input. Is there any way to accomplish this. I tried as follows:
process A{
input:
tuple val(some_name), path(some path) from some_channel
output:
tuple val(some val), path(some path) into A_some_ch
tuple va(some other va), path(some other path) into A_some_other_ch
val true into done_ch
script:
"""
some command out some files
some other command out some other files
"""
}
process B {
input:
val flag from done_ch
val (some val), path(some path) from A_some_ch
val (some other val), path(some other path) from A_some_other_ch
output:
some val and path
script:
"""
some commands out some files
"""
}

I got the idea using of flag from here: https://github.com/nextflow-io/patterns/blob/master/docs/mock-dependency.adoc
But if I use the flag, I got:
Operation 'set' does not allow more than one target name

Does it clarify the issue? Any Help?
Best, Zillur

Nextflow uses channels (asynchronous FIFO queues) to communicate between processes. The interaction between these processes, and ultimately the pipeline execution flow itself, is implicitly defined by these input and output declarations.
The issue with your code, I think, is the use of multiple input channels. Consider the following example that uses only a single queue channel. To avoid declaring the inputs again in process A's output declaration, we can output a key to join on as the first element in the output tuple. In the case below, this is simply the 'some_name' variable. We can then join up the two channels using the join operator as follows:
Channel
    .fromFilePairs( './data/*.txt', size:1 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .into { input_ch_A; input_ch_B }

process A {

    input:
    tuple val(some_name), path(some_path) from input_ch_A

    output:
    tuple val(some_name), val("procA_${some_name}"), path("procA_${some_name}") \
        into procA_results

    """
    touch "procA_${some_name}"
    """
}

process B {

    input:
    tuple val(some_name), path(some_path), val(some_other_val), path(some_other_path) \
        from input_ch_B.join(procA_results)

    """
    true
    """
}
```


text:
What is the exact definition of "breakend"?
I'm reading the Manta User Guide
→ it is a method to discover structural variants and indels from next generation sequencing data.
The word 'breakend' is used in several contexts
I have searched for a definition in duckduckgo, qwant and google but have not found a concise definition.
Thanks for you help

A structural variant is an inversion, deletion, duplication, or translocation. A "breakend" is one endpoint of a structural variant. This paper (Alkan et al.  Nature Review Genetics 2011) may help you understand: https://www.nature.com/articles/nrg2958
The Manta paper and documentation seem sparse on diagrams. DELLY, another structural variant caller, has diagrams in their publication. Their paper may be more of use: https://academic.oup.com/bioinformatics/article/28/18/i333/245403
Specific example: Say you're sequencing a chronic myeloid leukemia (CML) tumor. You will likely find reads that half map to the ABL1 gene on chr9 and half map to the BCR gene on chr22 (https://en.wikipedia.org/wiki/Philadelphia_chromosome). The genomic position of the location where chr9 turns to chr22 is called the breakpoint. The positions on chr9 and chr22 are the breakends.

Is it correct to say that the breakpoint is always on the sequence that I'm aligning to the reference, while the breakends always refer to the genomic positions in the reference genome?

I'd say breakpoint and breakends both refer to the genomic positions in the reference genome. A breakpoint typically contains two breakends, one for each genomic position.
The sequences you're aligning to the reference would be called "supporting reads". They could be split (half aligns to one position, half to another), or paired-end reads (one read maps to one position, its mate pair to another).

I think the term "breakend" is derived specifically from the VCF file spec:
https://samtools.github.io/hts-specs/VCFv4.2.pdf
Starting on page 12:
"5.4 Specifying complex rearrangements with breakends
An arbitrary rearrangement event can be summarized as a set of novel adjacencies. Each adjacency ties together 2 breakends. The two breakends at either end of a novel adjacency are called mates.
There is one line of VCF (i.e. one record) for each of the two breakends in a novel adjacency. A breakend record is identified with the tag “SVTYPE=BND” in the INFO field. The REF field of a breakend record indicates a base or
sequence s of bases beginning at position POS, as in all VCF records. The ALT field of a breakend record indicates a replacement for s. This “breakend replacement” has three parts:

The string t that replaces places s. The string t may be an extended version of s if some novel bases are inserted during the formation of the novel adjacency.

The position p of the mate breakend, indicated by a string of the form “chr:pos”. This is the location of the first mapped base in the piece being joined at this novel adjacency.

The direction that the joined sequence continues in, starting from p. This is indicated by the orientation of square brackets surrounding p. These 3 elements are combined in 4 possible ways to create the ALT. In each of the 4 cases, the assertion is that s is replaced with t, and then some piece starting at position p is joined to t. The cases are:

REF ALT Meaning
s t[p[ piece extending to the right of p is joined after t
s t]p] reverse comp piece extending left of p is joined after t
s ]p]t piece extending to the left of p is joined before t
s [p[t reverse comp piece extending right of p is joined before t
The example in Figure 1 shows a 3-break operation involving 6 breakends. It exemplifies all possible orientations of breakends in adjacencies. Notice how the ALT field expresses the orientation of the breakends.
#CHROM POS ID REF ALT QUAL FILTER INFO
2 321681 bnd W G G]17:198982] 6 PASS SVTYPE=BND
2 321682 bnd V T ]13:123456]T 6 PASS SVTYPE=BND
13 123456 bnd U C C[2:321682[ 6 PASS SVTYPE=BND
13 123457 bnd X A [17:198983[A 6 PASS SVTYPE=BND
17 198982 bnd Y A A]2:321681] 6 PASS SVTYPE=BND
17 198983 bnd Z C [13:123457[C 6 PASS SVTYPE=BND"


text:
GATK 4.2 using IndexFeatureFile with multiple files
I am trying to use GATK IndexFeatureFile to create ".vcf.gz.tbi" file. But I have 480 files.
For the pilot study I used:
!./gatk IndexFeatureFile -I HG00096.chrom22.bam.10061.vcf.gz -O HG00096.chrom22.bam.10061.vcf.gz.tbi

Doing this procedure for 480 files is not seem very practical. Is there any way to use all file as an input and it takes all files and gives ".tbi" files as an output?
I thought
!./gatk IndexFeatureFile -I *.vcf.gz -O *.vcf.gz.tbi

or maybe not using output argument (just input). I did not worked for me.

You could run loop over the files in your shell, running a separate command for each file. For example, in Bash (or another Bourne shell):
for file in *.vcf.gz; do ./gatk IndexFeatureFile -I ${file} -O ${file}.tbi; done

Is this approach also valid for merging vcf files? Because, normally I am using bcftools for merging. For example for 3 vcf files I use this: bcftools merge HG00096.chrom22.bam.10061.vcf.gz HG00097.chrom22.bam.10062.vcf.gz HG00101.chrom22.bam.10065.vcf.gz  -o merged.vcf

In the case of bcftools merge, I guess looping shouldn't be necessary since you can merge them with one command: bcftools merge -O z -o merged.vcf.gz HG*.vcf.gz


text:
nextflow: Filter outputs of a process
How could I filter outputs of a process in the input of the next process? Filtering works fine in channel, but if I try to filter outputs I got compilation error. I tried as follows:
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}
    process step1a {
    input:
    tuple val(sample_name), path(bfiles) from data
    output:
    tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)
    script:
    """
    plink --bfile "${sample_name}" --chr 1-22 --out "${prefix}.step1b" --make-bed
    """
    }
    process step1b {
    input:
    tuple val(bim), path(bims) from step1a_bim.filter{bim, files -> name =~/*bim/}
    tuple val(sample_name), path(bfiles) from step1a_out
    output:
    path("${prefix}.step1b.snplist.txt")
    tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
    script:
    """
    awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""

It gives me
Compilation error
-cause: Unexpected input: '{'
process step1b {

Step1a gives me chr1.step1b.bed|bim|fam, chr2.step1b.bed|bim|fam.......chr22.step1b.bed|bim|fam (22 plink files)
I need to access prefix of those outputs to run plink and other tools to generate chr1.step1c.bed|bim|fam ...chr22.step1c.bed|bim|fam (another 22 plink files). run awk to generate chr1.step1b.snplist.txt .... chr22.step1b.snplist.txt (22 txt files).
Any help to filter the inputs so that I could use only *bim files from step1a_bim? Also how to select prefix of files in the input of processes that are from step1a_out to run plink? Does it make sense now?
Best Reagrds
Zillur

In case of a compilation error it is more useful if you post your actual code, properly formatted, and your error message.

@Pallie, I edited my question above. Does it clarify my issues on outputs and inputs?

Instead of outputting all files and filtering for only bim, output only bims:
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)

Do
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.bim") into step1a_bim

So you can skip your complicated filter:
process step1b {
    input:
    tuple val(bim), path(bimfile) from step1a_bim
    output:
    path("${prefix}.step1b.snplist.txt")
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into
step1b_results
    script:
    """
    awk '{ if ((\\\\\$5=="T" && \\\\\$6=="A")||(\\\\\$5=="A" && \\\\\$6=="T")||(\\\\\$5=="C" && \\\\\$6=="G")||(\\\\\$5=="G" && \\\\\$6=="C")) print \\\\\$2, "ambig" ; else print \\\\\$2 ;}' "$bimfile" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"

"""

```

Thanks @Pallie . I tried your method. In that case, I got Process 'step1b' input file name collision -- There are multiple input files for each of the following file names: chr5.step1b.bim. Any suggestion?

The step1a_bim is a subset of step1a_out which produce bed|bim|fam files.

I could divide process step1b into two processes and could use step1a_out. But how could I use only *.bim files as input? Because step1a_out produce all *.bim, *.bed and *.fam files. For one process I need all three (bed|bim|fam) and for another I need only bim files. I am not sure how to subset the elements of a channel.


text:
Un-normalize DESeq2 counts
I obtained a matrix of RNA-seq count data that has been normalized by DESeq2's median of ratio method.
I know that DESeq2 wants to take in un-normalized counts, but I do not have access to those data.
How do I best proceed here if I want to perform DEG analysis using DESeq2?
I know I can always start from .fastq files, but that would be so much extra work.
I don't think I can un-normalize the data because I don't know the scaling factors.
Any suggestions?

You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?
It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.
I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)

Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:
estimateSizeFactorsForMatrix(round(ncounts))
 sample1  sample2  sample3  sample4  sample5  sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006

If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.
Below you can use the normalized counts and you can see the difference is not huge:
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)

Compare the new count table:
head(counts(dds,normalize=TRUE))
      sample1 sample2 sample3 sample4 sample5 sample6
gene1       0       2       0       1       0       0
gene2       9       8      16       6       1       2
gene3      23      40       8      38      11      41
gene4       0       0       0       6       8       0
gene5       3       2       0       5       1       2
gene6       0       1       3      19       4       9

And your input:
head(ncounts)
        sample1    sample2   sample3    sample4   sample5   sample6
gene1  0.000000  2.2658464  0.000000  0.9310803  0.000000  0.000000
gene2  8.693480  8.3081033 16.251453  5.5864817  1.121781  2.086619
gene3 23.472395 40.0299524  7.584012 38.1742919 11.217808 40.689076
gene4  0.000000  0.0000000  0.000000  5.5864817  7.852465  0.000000
gene5  3.477392  1.5105642  0.000000  4.6554014  1.121781  2.086619
gene6  0.000000  0.7552821  3.250291 18.6216058  4.487123  9.389787

We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))

Again, I would say I agree with @ATpoint comments that if you can, start from the fastq again.. don't work with data you have doubts about.

your answer are always  as always helpful..your suggestion or guidance on this https://bioinformatics.stackexchange.com/questions/15660/data-format-for-pathway-based-clustering-of-samples

You might be able to cheat by giving DESeq a custom size factor matrix where all the values are 1.  Then the software will feel like it's done the normalization step, but the 'normalized' values will be what you want.

I'll give this a try. Caveat is that I need to first round the normalized counts to integer, but I guess I'll have to live with that for now.

Nah, don't do that. This is most likely going to mess up things as the DESeq2 model runs on raw counts and everything (normalization, covariates) is done via offsets of the model. The DESeq2 author has been asked this a million times at Bioconductor support page, and if things were that easy then it would be in the vignette, at least as a quick'n'dirty workaround. Start from fastq, you are only adding uncertainty with these custom approaches that in the end cost you time and give inaccurate results. Learn from this, try to backup files and save your R environment incl. the raw data.


text:
RNASeq analysis using featureCount and EdgeR
I am using a pipeline (bam -> featurecount-> EdgeR) to do some RNASeq analysis of several groups and sub-groups. For example, I have the following dataset with two types (T1 and T2) and T1 has three groups G1, G2 and G3, T2 has G5 and G6 groups. I would like to find TPM numbers, fold-change and FDR for all the genes and Groups.
Type Group SRA_Run#
T1   G1.   <some SRA ids e.g. SRR123456, SRR112233>
T1   G2.   <some SRA>
T1.  G3.   <some SRA>
T2.  G5.   <some SRA>
T2.  G6.   <some SRA>

I have individual BAM files generated for each SRA Run. I also individual read counts from featureCount for each SRA Run. I was not sure if I need to merge individual read counts for each group or run featureCount using all bam files of a group at once. Does it make any difference in  the analysis that follows.
I also would like to know if I can run EdgeR at once for all the types and groups at once for doing differential gene expression for each group of Type 1 (T1) against every other group of Type 2 (T2). I am bit stuck at this point how to proceed. Any help would be appreciated.

It makes no difference if you process the BAM files one at a time with featureCounts or all together, except that it changes how you have to read the files into R.
You can supply edgeR with lists of contrasts to have it compute fold-changes and p-values for. Please have a look at the edgeR user guide for examples.

Thanks @Devon !! So, I can simply merge the count files, right? Also, I am taking fastq from different SRA runs, so they can be considered as replicates? I am asking because I have some dataset that contains 10 SRA runs for one sample and 4 SRA runs for another sample.

It depends on whether they're runs of the same sample or not. If they're from the same sample then they're not replicates (they can be merged together).

Thanks. Yes, they are from the same sample. In my original question, group G1 is the name of the sample and it has few SRA runs. Just one more question, will it be an issue if I use different number of SRA runs for two different groups? Also, one pipeline uses RUVr for normalizing the counts, so is it advised to run RUVr for merged featureCount results before running EdgeR? Please pardon my lack of knowledge about RNA-Seq analysis as I am bit new to this.

EdgeR and similar packages will handle the library size normalization. Nothing should be done on the individual runs before merging them.


text:
Telling grep to treat N as [ATCG]
Okay so I'm using grep to try and get a preview of some trimming operations that are not going as expected..
Lets say that my sequence in the FastQ file is: ATNGCNATCG
What I want to do is..
grep "ATCGCTATCG" my.fastq
..and match the sequence given above
Surely there is some way or some existing tool that I can use besides doing:
grep "[A|N][T|N][C|N][G|N][C|N]...etc."

Your fix is going to match NNNNN, ANNNC and etc. - are you sure this is what you want?

Technically, no it's not what I want, but as this is just trying to get a picture for what's going on under the hood I think it wont be a problem :)

Note that this will fail if the target sequence is split across multiple fastq lines. The fastq format does not require only one line of sequence per entry, you can have multiple ones.

If you want to stick to grep, use a scripting language such as Perl to generate the regex programmatically. For example:
perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' ATCGCTATCG

Prints:
[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]

You can use it in grep like so:
grep '[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]' <<< ATNGCNATCG

Prints:
ATNGCNATCG

If that works for you, you could make it into a little bash function that also runs the grep. Add these lines to your ~/.bashrc:
grepN(){
    seq="$1"
    file="$2"
    pattern=$(perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' "$seq")
    grep "$pattern" "$file"
}

You can now run
grepN ATCGCTATCG my.fastq

Of course, this is not a good idea since the sequence might be in different lines, but that's what you were doing originally.

Awesome, I'll try this out. It will at least save me the time typing out all those N's hah.

It does involve one extra step, but saves me tons of time and gets me the data I need! Thanks!

@RPINerd I edited Timur's very clever solution to add a shell function that you can use to do the whole thing in one step. I hope you don't mind, Timur, please let me know if you'd rather I posted a separate answer. It just seemed like a trivial extension of yours and you did all the work so...

@terdon Your edits and suggestions are always welcome! Thank you.