text:
Singularity And CNV Kit
I would like to try an use Singularity with CNV kit as the HPC I use seems to work better with Singularity:
Initial Attempt:
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "ls /"

Result:

bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var

Attempt 1:
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /proc/ && ls"

This results in a folder with a series of numbers but not the cnvkit py files.
Attempt 2:
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /singularity/ && ls"

This results in /bin/bash: line 1: cd: /singularity/: Not a directory
Attempt 3:
 singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /media/ && ls"

This results in nothing being printed to the screen. I tried this command three times to be sure.
How can I find the Singularity path I would need to use with singularity exec to run from the Singularity image.

I like this question a lot. Firstly I don't know singularity, but I know Docker and it will be similar. docker exec is very similar.
So for example the reasons are likely as follows:

media could be empty and hence no output.
the Python script may not be in /proc
singularity is likely to be a file rather than a directory.

One solution is,
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "ls -l /"

i.e. ask it which are the directories. Another solution is..
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "find / -name 'cnvkit*'"

This will locate cnvkit.py, note the use of ' and " for this sort of command.
Personally, I would dive into the container and issue the command of interest, manually. This way you get a feel for how it's constructed and don't need to keep issuing commands from outside the container. In Docker this is easy and it's impossible to keep a third party out of a container (even if the developer tries). In singularity it is likely similar.

From the notes I would simply access the singularity container directly.
Accessing a singularity container is different from Docker (I looked it up) ...
singularity shell ~/CNV_kit/cnvkit_latest.sif

Once inside you'll be able to look around, try out the find command, look in the bin etc... There will be no need to issue commands from outside the container to get it working.

From further comments,
A node on an HPC is a group of usually 8 to 32+ CPUs. What you have is a 'master' or 'head' node, where you log into and keep your data, and the stuff you do the work on is called 'nodes'. It has a more specific name, but it's not very cool, so we don't use it. The HPC will comprise a large number of nodes, e.g. 100+, but can be smaller particularly with the massive CPU machines around now. A node concept is a bit different on a NVidia machine if its part of a cluster.
The issue is someone could have called a script 'node', which would be confusing. It's complicated, because someone could have their own research machine which is attached and maintained on. the HPC, but only they can access and called it e.g. 'Joes_node', of 'Jills_node'.
Personally, I would sit down with your administrator and ask them to explain singularity and this the comments herein and this. They will know that at admin level (thats their job).
From my point of view, if I'd access to a cluster using singularity I could figure it out. I use GCP now, it's geeky but cool, so my information is out of date on non-Cloud clusters.

RamRS comments
What RamRS is saying and what I'm trying to assess whether you are within a singularity container. Again, I don't know the architecture.
Merging With Accepted Answer:
The main file in cnvkit is cnvkit.py in order to find that you need to look using singularity exec. The following command gives the path to the cnvkit.py file:
singularity exec cnvkit_latest.sif /bin/bash -c "cd /usr/local/bin/ && ls"

An example of how to use this path is below:
# Directory containing sample files:
directory="/path/to/bam/files/folder"

# Define path to singularity image:
cnvkit_container="/path/to/cnvkit_latest.sif"

cnvkit_path="/usr/local/bin/cnvkit.py"

singularity exec "$cnvkit_container" "$cnvkit_path" call "$directory/sample.cns" -o "$directory/sample.call.cns"

For some reason on my HPC it just says permission denied but it does not allow me to find the path?

I'd recommend looking under /bin/. Also, is there a Dockerfile like config that defines an entry point? If so, it should be straightforward to find out where things lie.

I can't access docker via my HPC. And under bin there is FileCheck-11         icuinfo     resolvectl  GET          id     rev  HEAD          infocmp     rgrep  POST          infotocap     rm but not anything that I can tell is directly related to the cnv kit py files?

@Indira Wait there's a misunderstanding singularity vs Docker. I would access the singularity container directly. The command to do so is below the ---- line in my post.

but wouldn't I need to use singularity exec to access the path to be able to use CNV kit?

Use singularity shell to find and test the executable, then use singularity exec to use it.

@Indira if you use the shell command it will give you access to the container - it operates just like an enclosed system (because it is). You will then be able to figure out how to run it via exec. If you try the shell command you will see after e.g. pwd.

@RamRS beat me to it and yes 100% agree.

@M__ I was unaware that singularity shell is a thing, thank you so much!

I did singularity shell cnvkit_latest.sif and then find /path/to/cnvkit_latest.sif -type f -name "*.py" but it doesn't give any output?

Run find / -name "c*.py". You're within the singularity container's shell, so searching with the path to the container does not make sense.

Even with singularity shell cnvkit_latest.sif when I run find / -name "c*.py". it doesn't work - it searches not just in the singularity but through all the nodes that other users use. I did get one path that seemed like what I would be looking for '/etc/ssl/private' and it said Permission denied so I tried singularity exec cnvkit_latest.sif /bin/bash -c "cd /etc/ssl/private/ && ls" and it still said Permission denied

Are you sure you're running the command from within the singularity shell? Compare the output to singularity exec file.sif /bin/bash -c "ls /" run from your HPC shell and a simple ls / from within singularity shell. The output should be identical

I am sure I am running it in the shell

@Indira what is the output of pwd after the shell command?

I told that the shell command with the find searches through other people nodes and pulls up their files with a similar python files such as under Amy_node and Barbara_node, but I only want it under my node.

I see ... I know what you are doing via singularity  but I lack specific knowledge

By "node" do you mean compute node?

Yes I am one an hpc so it does mean a computing node.

Thank you know more than me about how to use this site.

I'll convert the comments to chat or a wiki shortly. The chat site is here: https://chat.stackexchange.com/rooms/58858/bioinformatics if you call a user e.g. @M__ that will trigger a notification to that member.

@M__ The nomenclature is now head and worker nodes. Even GitHub has now renamed the principal code branch to main, so the nomenclature you're censoring is outdated.

Is it possible to merge @M__ answer and my answer as it is on the same topic?

Thanks to OP for sharing that their image is from https://hub.docker.com/r/etal/cnvkit
I downloaded this and it seems pretty simple.
My tests:
singularity pull docker://etal/cnvkit:latest
$ singularity shell ./cnvkit_latest.sif
INFO:    Converting SIF file to temporary sandbox...
bash: warning: setlocale: LC_ALL: cannot change locale (en_US.UTF-8)
Singularity> ls /
bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var
Singularity> ls /singularity
/singularity
Singularity> cnvkit.py
usage: cnvkit.py [-h]
                 {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                 ...

CNVkit, a command-line toolkit for copy number analysis.

positional arguments:
  {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                        Sub-commands (use with -h for more info)
    batch               Run the complete CNVkit pipeline on one or more BAM files.
    target              Transform bait intervals into targets more suitable for CNVkit.
    access              List the locations of accessible sequence regions in a FASTA file.
    antitarget          Derive off-target ("antitarget") bins from target regions.
    autobin             Quickly calculate reasonable bin sizes from BAM read counts.
    coverage            Calculate coverage in the given regions from BAM read depths.
    reference           Compile a coverage reference from the given files (normal samples).
    fix                 Combine target and antitarget coverages and correct for biases. Adjust raw coverage data according to the given reference, correct potential biases and re-center.
    segment             Infer copy number segments from the given coverage table.
    call                Call copy number variants from segmented log2 ratios.
    diagram             Draw copy number (log2 coverages, segments) on chromosomes as a diagram. If both the raw probes and segments are given, show them side-by-side on each chromosome (segments on the left side, probes on the right
                        side).
    scatter             Plot probe log2 coverages and segmentation calls together.
    heatmap             Plot copy number for multiple samples as a heatmap.
    breaks              List the targeted genes in which a copy number breakpoint occurs.
    genemetrics         Identify targeted genes with copy number gain or loss.
    sex                 Guess samples' sex from the relative coverage of chromosomes X and Y.
    metrics             Compute coverage deviations and other metrics for self-evaluation.
    segmetrics          Compute segment-level metrics from bin-level log2 ratios.
    bintest             Test for single-bin copy number alterations.
    import-picard       Convert Picard CalculateHsMetrics tabular output to CNVkit .cnn files. The input file is generated by the PER_TARGET_COVERAGE option in the CalculateHsMetrics script in Picard tools. If 'antitarget' is in the
                        input filename, the generated output filename will have the suffix '.antitargetcoverage.cnn', otherwise '.targetcoverage.cnn'.
    import-seg          Convert a SEG file to CNVkit .cns files.
    import-theta        Convert THetA output to a BED-like, CNVkit-like tabular format. Equivalently, use the THetA results file to convert CNVkit .cns segments to integer copy number calls.
    import-rna          Convert a cohort of per-gene log2 ratios to CNVkit .cnr format.
    export              Convert CNVkit output files to another format.
    version             Display this program's version.

options:
  -h, --help            show this help message and exit

See the online manual for details: https://cnvkit.readthedocs.io
Singularity> which cnvkit.py
/usr/local/bin/cnvkit.py

As you can see, since cnvkit.py is in /usr/local/bin/, simply typing cnvkit.py within the shell works. Let's try an exec with no frills:
$ singularity exec cnvkit_latest.sif cnvkit.py
INFO:    Converting SIF file to temporary sandbox...
usage: cnvkit.py [-h]
                 {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                 ...

CNVkit, a command-line toolkit for copy number analysis.

positional arguments:
  {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                        Sub-commands (use with -h for more info)
    batch               Run the complete CNVkit pipeline on one or more BAM files.
    target              Transform bait intervals into targets more suitable for CNVkit.
    access              List the locations of accessible sequence regions in a FASTA file.
    antitarget          Derive off-target ("antitarget") bins from target regions.
    autobin             Quickly calculate reasonable bin sizes from BAM read counts.
    coverage            Calculate coverage in the given regions from BAM read depths.
    reference           Compile a coverage reference from the given files (normal samples).
    fix                 Combine target and antitarget coverages and correct for biases. Adjust raw coverage data according to the given reference, correct potential biases and re-center.
    segment             Infer copy number segments from the given coverage table.
    call                Call copy number variants from segmented log2 ratios.
    diagram             Draw copy number (log2 coverages, segments) on chromosomes as a diagram. If both the raw probes and segments are given, show them side-by-side on each chromosome (segments on the left side, probes on the right
                        side).
    scatter             Plot probe log2 coverages and segmentation calls together.
    heatmap             Plot copy number for multiple samples as a heatmap.
    breaks              List the targeted genes in which a copy number breakpoint occurs.
    genemetrics         Identify targeted genes with copy number gain or loss.
    sex                 Guess samples' sex from the relative coverage of chromosomes X and Y.
    metrics             Compute coverage deviations and other metrics for self-evaluation.
    segmetrics          Compute segment-level metrics from bin-level log2 ratios.
    bintest             Test for single-bin copy number alterations.
    import-picard       Convert Picard CalculateHsMetrics tabular output to CNVkit .cnn files. The input file is generated by the PER_TARGET_COVERAGE option in the CalculateHsMetrics script in Picard tools. If 'antitarget' is in the
                        input filename, the generated output filename will have the suffix '.antitargetcoverage.cnn', otherwise '.targetcoverage.cnn'.
    import-seg          Convert a SEG file to CNVkit .cns files.
    import-theta        Convert THetA output to a BED-like, CNVkit-like tabular format. Equivalently, use the THetA results file to convert CNVkit .cns segments to integer copy number calls.
    import-rna          Convert a cohort of per-gene log2 ratios to CNVkit .cnr format.
    export              Convert CNVkit output files to another format.
    version             Display this program's version.

options:
  -h, --help            show this help message and exit

See the online manual for details: https://cnvkit.readthedocs.io
INFO:    Cleaning up image...

The answer seems to be a lot easier than having to find cnvkit.py, as it is already in the PATH.

This did not work for me as it seems to work differently on the HPC I am on.

That does not make much sense. There could be something unusual with your HPC. Perhaps you should contact your sysadmin.

I have done so and I simply have to navigate around these issues using the path inside the singularity image.


text:
Calculation of Site Frequency Spectrum (SFS) using STR data
Hi everyone I had a question on if it is possible to calculate (Site Frequency Spectrum)SFS using STR data. Please mention if you know any programm that can be used for that!
I have tried Arlequin but it doesn't offer this option for STR. Personally I haven't find anything that can do this type of calculation with STR!
Thank you!

@M__ "I have tried Arlequin but it doesn't offer this option for STR"

There are two things. First STS - which is a complex algorithm, secondly is the data STR. If Arlequin has the maximum likelihood estimate for STS (which the OP  implying it does), it's easy. This is because Arlequin will process any data format available, STR is simply a data format, not an algorithm. It's about looking at what format would be most suitable and converting the data accordingly to a format Arlequin uses. The Arlequin data types are not rigid like phylogenetics. That is a fairly easy question. If Arlequin doesn't have STS (the algorithm), finding an STS algorithms is not trivial.

The parsing operation might be a little complicated - but it's not hard and someone has probably written a parsing script, not necessarily STR -> Arlequin, but STR to an intermediate format that will pass into the official Arlequin data parser.

source: Stack Exchange

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

sort: votes


text:
NA values from a URL import
Trying to load the XML data from this url
URL: https://www.senate.gov/general/contact_information/senators_cfm.xml.
# Load required library
library(RCurl)

# Define the URL
url <- "https://www.senate.gov/general/contact_information/senators_cfm.xml"

# Retrieve XML content from the URL
xml_content <- getURL(url)

# Load XML library
library(XML)

# Parse the XML content
xml_parsed <- xmlParse(xml_content)

# Extract data from XML and convert it to a data frame
data <- xmlToDataFrame(xml_parsed)

# Display the dimensions of the data
dim(data)

But my last row is all NA's. I know the dim is suppose to be 100 due to the number of senators per state but it's currently 101 due to that empty row.
> data[1:5]
            member_full    last_name     first_name party state
1        Baldwin (D-WI)      Baldwin          Tammy     D    WI
2       Barrasso (R-WY)     Barrasso           John     R    WY
3         Bennet (D-CO)       Bennet     Michael F.     D    CO
4      Blackburn (R-TN)    Blackburn         Marsha     R    TN
5     Blumenthal (D-CT)   Blumenthal        Richard     D    CT
6         Booker (D-NJ)       Booker        Cory A.     D    NJ
7        Boozman (R-AR)      Boozman           John     R    AR
8          Braun (R-IN)        Braun           Mike     R    IN
9          Britt (R-AL)        Britt     Katie Boyd     R    AL
10         Brown (D-OH)        Brown        Sherrod     D    OH
11          Budd (R-NC)         Budd            Ted     R    NC
12      Cantwell (D-WA)     Cantwell          Maria     D    WA
13        Capito (R-WV)       Capito  Shelley Moore     R    WV
14        Cardin (D-MD)       Cardin    Benjamin L.     D    MD
15        Carper (D-DE)       Carper      Thomas R.     D    DE
16         Casey (D-PA)        Casey Robert P., Jr.     D    PA
17       Cassidy (R-LA)      Cassidy           Bill     R    LA
18       Collins (R-ME)      Collins       Susan M.     R    ME
19         Coons (D-DE)        Coons Christopher A.     D    DE
20        Cornyn (R-TX)       Cornyn           John     R    TX
21  Cortez Masto (D-NV) Cortez Masto      Catherine     D    NV
22        Cotton (R-AR)       Cotton            Tom     R    AR
23        Cramer (R-ND)       Cramer          Kevin     R    ND
24         Crapo (R-ID)        Crapo           Mike     R    ID
25          Cruz (R-TX)         Cruz            Ted     R    TX
26        Daines (R-MT)       Daines          Steve     R    MT
27     Duckworth (D-IL)    Duckworth          Tammy     D    IL
28        Durbin (D-IL)       Durbin     Richard J.     D    IL
29         Ernst (R-IA)        Ernst           Joni     R    IA
30     Feinstein (D-CA)    Feinstein         Dianne     D    CA
31     Fetterman (D-PA)    Fetterman           John     D    PA
32       Fischer (R-NE)      Fischer            Deb     R    NE
33    Gillibrand (D-NY)   Gillibrand     Kirsten E.     D    NY
34        Graham (R-SC)       Graham        Lindsey     R    SC
35      Grassley (R-IA)     Grassley          Chuck     R    IA
36       Hagerty (R-TN)      Hagerty           Bill     R    TN
37        Hassan (D-NH)       Hassan  Margaret Wood     D    NH
38        Hawley (R-MO)       Hawley           Josh     R    MO
39      Heinrich (D-NM)     Heinrich         Martin     D    NM
40  Hickenlooper (D-CO) Hickenlooper        John W.     D    CO
41        Hirono (D-HI)       Hirono       Mazie K.     D    HI
42        Hoeven (R-ND)       Hoeven           John     R    ND
43    Hyde-Smith (R-MS)   Hyde-Smith          Cindy     R    MS
44       Johnson (R-WI)      Johnson            Ron     R    WI
45         Kaine (D-VA)        Kaine            Tim     D    VA
46         Kelly (D-AZ)        Kelly           Mark     D    AZ
47       Kennedy (R-LA)      Kennedy           John     R    LA
48          King (I-ME)         King  Angus S., Jr.     I    ME
49     Klobuchar (D-MN)    Klobuchar            Amy     D    MN
50      Lankford (R-OK)     Lankford          James     R    OK
51           Lee (R-UT)          Lee           Mike     R    UT
52        LujÃ¡n (D-NM)       LujÃ¡n        Ben Ray     D    NM
53        Lummis (R-WY)       Lummis     Cynthia M.     R    WY
54       Manchin (D-WV)      Manchin       Joe, III     D    WV
55        Markey (D-MA)       Markey      Edward J.     D    MA
56      Marshall (R-KS)     Marshall          Roger     R    KS
57     McConnell (R-KY)    McConnell          Mitch     R    KY
58      Menendez (D-NJ)     Menendez         Robert     D    NJ
59       Merkley (D-OR)      Merkley           Jeff     D    OR
60         Moran (R-KS)        Moran          Jerry     R    KS
61        Mullin (R-OK)       Mullin      Markwayne     R    OK
62     Murkowski (R-AK)    Murkowski           Lisa     R    AK
63        Murphy (D-CT)       Murphy    Christopher     D    CT
64        Murray (D-WA)       Murray          Patty     D    WA
65        Ossoff (D-GA)       Ossoff            Jon     D    GA
66       Padilla (D-CA)      Padilla           Alex     D    CA
67          Paul (R-KY)         Paul           Rand     R    KY
68        Peters (D-MI)       Peters        Gary C.     D    MI
69          Reed (D-RI)         Reed           Jack     D    RI
70      Ricketts (R-NE)     Ricketts           Pete     R    NE
71         Risch (R-ID)        Risch       James E.     R    ID
72        Romney (R-UT)       Romney           Mitt     R    UT
73         Rosen (D-NV)        Rosen          Jacky     D    NV
74        Rounds (R-SD)       Rounds           Mike     R    SD
75         Rubio (R-FL)        Rubio          Marco     R    FL
76       Sanders (I-VT)      Sanders        Bernard     I    VT
77        Schatz (D-HI)       Schatz          Brian     D    HI
78       Schmitt (R-MO)      Schmitt           Eric     R    MO
79       Schumer (D-NY)      Schumer     Charles E.     D    NY
80         Scott (R-FL)        Scott           Rick     R    FL
81         Scott (R-SC)        Scott            Tim     R    SC
82       Shaheen (D-NH)      Shaheen         Jeanne     D    NH
83        Sinema (I-AZ)       Sinema        Kyrsten     I    AZ
84         Smith (D-MN)        Smith           Tina     D    MN
85      Stabenow (D-MI)     Stabenow         Debbie     D    MI
86      Sullivan (R-AK)     Sullivan            Dan     R    AK
87        Tester (D-MT)       Tester            Jon     D    MT
88         Thune (R-SD)        Thune           John     R    SD
89        Tillis (R-NC)       Tillis           Thom     R    NC
90    Tuberville (R-AL)   Tuberville          Tommy     R    AL
91    Van Hollen (D-MD)   Van Hollen          Chris     D    MD
92         Vance (R-OH)        Vance           J.D.     R    OH
93        Warner (D-VA)       Warner        Mark R.     D    VA
94       Warnock (D-GA)      Warnock     Raphael G.     D    GA
95        Warren (D-MA)       Warren      Elizabeth     D    MA
96         Welch (D-VT)        Welch          Peter     D    VT
97    Whitehouse (D-RI)   Whitehouse        Sheldon     D    RI
98        Wicker (R-MS)       Wicker       Roger F.     R    MS
99         Wyden (D-OR)        Wyden            Ron     D    OR
100        Young (R-IN)        Young           Todd     R    IN
101                <NA>         <NA>           <NA>  <NA>  <NA>

This is political data and might be more appropriate e.g. for https://stat.stackexchange.com.

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


text:
How Can You Reformat an OTU Table from qiime2 pathway to Include Sequence IDs?
I am trying to run Tax4Fun in MicrobiomeAnalyst, but am running into trouble with my table formatting. From the qiime2 pathway, I have an OTU, metadata and taxonomic table formatted like so:

In order to use Tax4Fun, I understand that I will need to reformat my OTU table to include sequence IDs so it looks like this example:

However, I cannot figure out how to merge my sequence files from the qiime2 pathway with the OTU table so it looks like that. Any suggestions would be much appreciated!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.


text:
Making a list of lists in R
I have a data frame columns are pathways and genes are rows
Like
head(a[1:2,1:2])
  Release.of.Cancer.Cell.Antigens Cancer.Antigen.Presentation
1                          ARID1A                        ARG2
2                             ATM                        ATF3

str(a[1:2,1:2])
'data.frame':   2 obs. of  2 variables:
 $ Release.of.Cancer.Cell.Antigens: chr  "ARID1A" "ATM"
 $ Cancer.Antigen.Presentation    : chr  "ARG2" "ATF3"

I want to make a list of pathways containing their genes, something like
str(h.sets)
List of 50
 $ HALLMARK_ADIPOGENESIS                     : chr [1:208] "ABCA1" "ABCB8" "ACAA2" "ACADL" ...
 $ HALLMARK_ALLOGRAFT_REJECTION              : chr [1:332] "AARS1" "ABCE1" "ABI1" "ACHE" ...

I done so
h.names <- unique(colnames(a))
h.sets <- vector("list",length=length(h.names))
names(h.sets) <- h.names

But here I get error here
for (i in names(h.sets)) {
    h.sets[[i]] <- pull(a[colnames(a)==i])
}

Error in `pull()`:
! Can't extract columns past the end.
ℹ Location 1 doesn't exist.
ℹ There are only 0 columns.
Run `rlang::last_trace()` to see where the error occurred.

Any help please

Data frame is a type of list, so for most purposes you already have what you're asking. You can convert it into an ordinary list with as.list(a). If that doesn't help, you should tell more about what you're trying to achieve.

Why are you using pull? Why not a simple h.sets[[i]] <- a[i]?


text:
Which type of variant caller should I use in a WES normal cell line sample?
I have whole-exome sequencing data of an immortalised non-tumor (normal) cell line that I wish to assess for the presence/absence of APC/Wnt mutations. This is to double check that the cell line is sufficiently "clean", without mutations in key regions, for further experimental analysis.
I understand that without a matched normal sample, I won't be able to distinguish if a variant is a somatic variant or a rare inherited germline mutation. As I just want to check for mutations in those regions, it is not important to distinguish between these types.
Question: Should I use the GATK haplotypecaller or a somatic variant caller such as Mutect2 without a matched-normal sample?

I don't understand. If this is not a tumor sample, why would you want to run a somatic variant caller? Do you want to run this in a super sensitive way to detect everything that might conceivably be a variant?


text:
Prokka + table2asn (NCBI submission)
I have a couple of MAGs of interest that I am trying to submit to the NCBI (with annotations).
The process of submitting the annotations is rather challenging to me. I have the output of Prokka, which includes multiple files, but the one I think I need is the .gff (which is GFF3 format).
The NCBI website descibes how to obtain the required file format (.sqn) from .gff files (https://www.ncbi.nlm.nih.gov/genbank/genomes_gff/). I downloaded the required file (table2asn), and ran the following code:
"table2asn" -M n -Z -J -c w -t "template.sbt" -V b -j "[organism=Organism] [strain=strain] [gcode=11]" -i "fasta.fna" -f "annotations.gff" -o output.sqn

However, I got an error:
GFF annotation requires locus tags, which are missing from one or more genes, so the command line argument -locus-tag-prefix is needed. Which was a bit strange, since I could see that the original file had a locus_tag entry. The format liiked kind of like this:
S6C7    Prodigal:002006 CDS 244817  246085  .   -   0   ID=AJNLIECJ_06723;inference=ab initio prediction:Prodigal:002006,similar to AA sequence:ISfinder:ISNpu7;locus_tag=AJNLIECJ_06723;product=IS701 family transposase ISNpu7
S6C7    Prodigal:002006 CDS 246142  246702  .   -   0   ID=AJNLIECJ_06724;inference=ab initio prediction:Prodigal:002006;locus_tag=AJNLIECJ_06724;product=hypothetical protein

Anyway, I added -locus-tag-prefix AJNLIECJ to my code, which at least got me an output, but now the output is a touch weird:
 gene            complement(246142..246702)
                 /locus_tag="AJNLIECJ_006724"
 CDS             complement(246142..246702)
                 /locus_tag="AJNLIECJ_006724"
                 /inference="ab initio prediction:Prodigal:002006"
                 /codon_start=1
                 /transl_table=11
                 /product="hypothetical protein"
                 /protein_id="AJNLIECJ:AJNLIECJ_06724"

So basically, the /protein_id is repeated, and sometimes the numbers between locus and protein ID don't even match (some of the protein ids 'go ahead' for some reason).
My thoughts:

I have no idea how to fix this
Is it even necessary to fix it? Since all NCBI protein IDs tend to be in a consistent format anyway, so no matter what I add, it will get updated.
Should I remove the /protein_id field? How?

Thank you for your time!

Hi @Laura, contact the NCBI Helpdesk. They will respond with the exact answer within 1 day and can be within a few hours. You can copy paste, or just cite this question on the request form.


text:
Nucleotide alignments from mmseqs "tblastn"
I am searching a protein sequences against a nucleotide database ("tblastn") using mmseqs and need nucleotide sequence alignments as an output.
I start by the following (pg 23 of https://mmseqs.com/latest/userguide.pdf)
mmseqs search ecoli_proteins ecoli_genome alnDB tmp

and then get alignments as follows:
mmseqs convertalis ecoli_proteins ecoli_genome alnDB alnDB.html --format-mode 3

However, this gives me protein sequence alignments, when I really want the nucleotide base alignments (i.e. would be aligned codon-by-codon).
What is the proper way to get this output from mmseqs?


text:
Integrate bulk RNA and ATAC-seq genes
I have up and down-regulated genes from bulk RNA result from DeSeq2. Genes that are differential in chromatin accessibility between two conditions from DiffBind. I would like to find the genes overlap between bulk RNA and ATAC-seq assay as these papers but don't know the code to get it. Would you please have a suggestion? Thank you so much  .

paper 1
paper 2

in silico DEG details
I used intersect function in R which can do the same job. Which I am not sure is the experiment try to find transcription factors that cause the difference in phenotype between 2 conditions. DEG in bulk-RNA from DeSeq2 (log2 fold change >2.5 or log2 fold change < -2.5) intersect with genes differently in chromatin accessibility from DiffBind (log2 fold change >2.5 or log2 fold change < -2.5). I narrowed down to 142 genes.
For atac-seq result, after applying filter log2fold > 2.5 or log2fold < -2.5 and p value < 0.00001, I still got around 12k genes that are different between 2 conditions. Meanwhile bulk-RNA after apply filter has only around 300 genes. The number of genes in human and pig are similar so I am not sure about the number of gene I have in atac-seq
What do you think about this approach because I don't use MEME or TOMTOM?

I share with you what the author did.

" In order to screen new transcription factors by combining with star transcription factors. First, we used MEME software to predict the motif of Duroc pigs' peaks. Second, we compared the transcription factor database with TOMTOM to predict all transcription factors that may bind to the chromatin open regions . Third, we found all possible binding sites of MYOD, MYOG, MYF5, MYF6 and MEF2 families and did a union. A total of 99 sites were found.
"Fourthly, we counted other transcription factors that may bind to these 99 sites and sorted them by the number of binding sites. All of the above statistics were manually counted."

Will read more to give you a good reply. Thank you so much

You requested I make a response. Okay, I'm reluctant to do this firstly, all my  eukaryotes genomics is out of date except for immunology stuff and secondly, this is a hardcore eukaryote gene expression question. I am really not sure about the pig-human side of this experiment.
What you have is an extremely good data set. What you appear to be doing is assessing the complexity of chromatin based gene expression, using very hardcore data. What we don't know is the variation e.g. in experiments leading to the RNA-seq. I would guess the RNA-seq is varying and being compared against atac-seq.
I have not read the papers, read the abstracts you posted.
MEME and TOMTOM don't seem to be impressive in comparison to an absolute biological experiment RNA-seq vs. atac-seq. Predictions are never perfect (the very, very best predictions are 80%) and you have an absolute data set - why even be concerned about predictions based analysis?
I do know DEG and the approach you are using is very conservative. It's good to visualise this on a volcano plot. Personally I just use a probabilistic threshold and then assess the log2 for each gene group. You are using both -/+2.5 log2 and <1E-6 . Personally I would be concerned with the stringency and also look at results with e.g. just <1E-6 or <1E-6 and -/+ 0.5log2. Keep in mind I don't do eukaryotes.
The threshold criteria will be the number 1 issue that affects your results and that is the most important information I can give. So do your approach and also try a more relaxed threshold. Personally I wouldn't relax the p<1E-6, but I would relax the log2. The reasons for this are a bit complicated, but they are statistically sound.
Your intersect approach seems fine, viz.
intersect(object1, object2)

@haci's method is cool, the above does the job too. The question is the genes vs thresholds and the experimental set ups and then the comparison between experimental setups, i.e. which genes are switching under chromatin control.
What I am not clear about is the experimental design and how you are performing DSeq2,

I assume you are using an internal control for RNAseq and a comparable one for atac-seq. I assume you are not comparing atac-seq verse RNAseq in the DEG.
The alternative is atac-seq vs. RNAseq DEG and the comparison is pig versus human. I honestly dunno.

Anway the caveat I would add is an intersect SHOULD NOT ignore whether it is +log2 or -log2 ... you need to do four intersects

+ve verse +ve
-v verse -ve
+ve verse -ve
-ve verse +ve

Again its four separate intersects.
Good luck looks a really interesting experiment and analytics look simple but powerful in context.

Thanks @M__  so much for a very detailed reply! Each RNA-seq and atac-seq has control vs diseased. I can do a volcano plot. So I compared RNA-seq control vs RNA-seq diseased. ATAC seq control vs ATAC seq diseased. Then I try to figure out which genes/transcription factors may be the cause of the difference in the phenotype between two conditions. I see papers do scRNA-seq integrative ATAC-seq but not many papers do bulk RNA-seq integrative ATAC-seq. Maybe I should intersect up-regulated genes vs open chromatin region genes and down-regulated genes vs close chromatin region genes.

Yes, I tried to do it at the beginning but I need two more point to allow upvote. Sorry for that! When I have enough reputation, I will upvote all your answers right away!

I am allowed to upvote now .

I've borrowed one of your smilies, thank you @chris

You are welcome .

If I am getting your question right, all you need is https://bioinfogp.cnb.csic.es/tools/venny/. I assume that your two lists corresponding to RNA-seq and ATAC-seq use the same set of gene identifiers (gene symbol, Entrez ID, Ensembl ID etc).
The Results section is populated with entries that belongs to the corresponding "region" upon clicking on the regions of the Venn diagram.
Personally I have never used it but my wet-lab colleagues use it all the time.


text:
ProteinPromt installation
i dont know how to build/pull this container (https://codebase.helmholtz.cloud/proteinprompt/ProteinPrompt) can you advice me how to install it?
EDIT:
My soft - Ubuntu 20.04 lts
Thanks to @M__ I manage to build image (almost). But I have an issue :
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker build -t dockerfile .
[+] Building 1.0s (15/15) FINISHED
 => [internal] load build definition from Dockerfile                       0.0s
 => => transferring dockerfile: 1.19kB                                     0.0s
 => [internal] load .dockerignore                                          0.0s
 => => transferring context: 2B                                            0.0s
 => [internal] load metadata for docker.io/continuumio/miniconda3:4.10.3-  0.9s
 => [ 1/10] FROM docker.io/continuumio/miniconda3:4.10.3-alpine@sha256:e0  0.0s
 => [internal] load build context                                          0.0s
 => => transferring context: 650B                                          0.0s
 => CACHED [ 2/10] RUN rm /bin/sh && ln -s /bin/bash /bin/sh               0.0s
 => CACHED [ 3/10] RUN conda update -n base conda                          0.0s
 => CACHED [ 4/10] COPY ./data/proteinPrompt.yml /                         0.0s
 => CACHED [ 5/10] RUN conda env update --name root -f /proteinPrompt.yml  0.0s
 => CACHED [ 6/10] RUN mkdir /proteinPrompt/                               0.0s
 => CACHED [ 7/10] COPY ./proteinPrompt/*py /proteinPrompt/                0.0s
 => CACHED [ 8/10] RUN mkdir /model/                                       0.0s
 => CACHED [ 9/10] RUN wget https://zenodo.org/record/5708224/files/rf_hu  0.0s
 => CACHED [10/10] COPY ./data/all.proteins.HUMAN.autocorrelation.pbz2 /m  0.0s
 => exporting to image                                                     0.0s
 => => exporting layers                                                    0.0s
 => => writing image sha256:ae93ff408fb6468280a5f0344bb2ac5ac811df2e5d32a  0.0s
 => => naming to docker.io/library/dockerfile                              0.0s
WARNING: buildx: failed to read current commit information with git rev-parse --is-inside-work-tree

and image is named 'dockerfile'? and it do not work
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker images
REPOSITORY    TAG           IMAGE ID       CREATED          SIZE
dockerfile    latest        ae93ff408fb6   11 minutes ago   2.02GB
alphafold     latest        c41a4b5e9d48   13 days ago      10GB
nvidia/cuda   11.0.3-base   4f3e95122fab   23 months ago    122MB

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker run  -v $(pwd)/data:/data proteinprompt search -f data/input.fa -o data/output.csv
Unable to find image 'proteinprompt:latest' locally
docker: Error response from daemon: pull access denied for proteinprompt, repository does not exist or may require 'docker login': denied: requested access to the resource is denied.
See 'docker run --help'.

EDIT 2:
I try to edit dockerfile (in several places)

FROM continuumio/miniconda3:4.10.3-alpine
MAINTAINER Sebastian Canzler <sebastian.canzler@ufz.de>
LABEL authors="sebastian.canzler@ufz.de" \
    description="Docker image of ProteinPrompt - a tool to predict protein-protein interactions based on amino acid sequences."
RUN apt-get update && apt-get install -y git #<-----HERE

#Try This to get a bash in order to run 'source activate modfinder'
RUN rm /bin/sh && ln -s /bin/bash /bin/sh

RUN conda update -n base conda
COPY ./data/proteinPrompt.yml /
RUN conda env update --name root -f /proteinPrompt.yml && \
    rm -rf /opt/conda/pkgs/*
ENV PATH /opt/conda/envs/proteinPrompt/bin:$PATH
# copy the proteinPrompt files
RUN mkdir /proteinPrompt/
COPY ./proteinPrompt/*py /proteinPrompt/
ENV PATH /proteinPrompt:$PATH

# get the random forest model from zenodo
RUN mkdir /model/
#COPY ./data/rf_human_750_inv.rev.pkl /model/random_forest.pkl
RUN wget https://zenodo.org/record/5708224/files/rf_human_750_inv.rev.pkl?download=1 -O /model/random_forest.pkl
COPY ./data/all.proteins.HUMAN.autocorrelation.pbz2 /model/all.proteins.HUMAN.autocorrelation.pbz2

# make proteinPrompt executable
ENTRYPOINT ["python", "/proteinPrompt/ProteinPrompt.py"]
#CMD ["-h"]

but there is a bash problem:
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker build -t dockerfile .
[+] Building 1.6s (10/15)
 => [internal] load .dockerignore                                          0.0s
 => => transferring context: 2B                                            0.0s
 => [internal] load build definition from Dockerfile                       0.0s
 => => transferring dockerfile: 1.23kB                                     0.0s
 => [internal] load metadata for docker.io/continuumio/miniconda3:4.10.3-  1.0s
 => [ 1/11] FROM docker.io/continuumio/miniconda3:4.10.3-alpine@sha256:e0  0.0s
 => [internal] load build context                                          0.0s
 => => transferring context: 650B                                          0.0s
 => CACHED [ 2/11] RUN rm /bin/sh && ln -s /bin/bash /bin/sh               0.0s
 => CACHED [ 3/11] RUN conda update -n base conda                          0.0s
 => CACHED [ 4/11] COPY ./data/proteinPrompt.yml /                         0.0s
 => CACHED [ 5/11] RUN conda env update --name root -f /proteinPrompt.yml  0.0s
 => ERROR [ 6/11] RUN apt-get update && apt-get install -y git             0.5s
------
 > [ 6/11] RUN apt-get update && apt-get install -y git:
#0 0.537 /bin/ash: apt-get: not found
------
WARNING: buildx: failed to read current commit information with git rev-parse --is-inside-work-tree
Dockerfile:15
--------------------
  13 |         rm -rf /opt/conda/pkgs/*
  14 |     ENV PATH /opt/conda/envs/proteinPrompt/bin:$PATH
  15 | >>> RUN apt-get update && apt-get install -y git
  16 |
  17 |     # copy the proteinPrompt files
--------------------
ERROR: failed to solve: process "/bin/ash -eo pipefail -c apt-get update && apt-get install -y git" did not complete successfully: exit code: 127

EDIT 3:
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker pull registry.hzdr.de/proteinprompt/proteinprompt:latest
latest: Pulling from proteinprompt/proteinprompt
5843afab3874: Already exists
72c585c386a3: Already exists
87d798ddb41c: Already exists
316f048b0253: Pull complete
ae11eb6b6512: Pull complete
d3dc3767a31c: Pull complete
ea3dbaf26ce1: Pull complete
061c354f8f77: Pull complete
60eb5484e32a: Pull complete
0625cf0fe317: Pull complete
442818fdc663: Pull complete
203550cafc1f: Pull complete
Digest: sha256:01f0251f34313437e5fcc2dbb26d1ec1103d50aef6335b2261a2048be8ed0813
Status: Downloaded newer image for registry.hzdr.de/proteinprompt/proteinprompt:latest
registry.hzdr.de/proteinprompt/proteinprompt:latest

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker images
REPOSITORY                                     TAG           IMAGE ID       CREATED         SIZE
alphafold                                      latest        c41a4b5e9d48   13 days ago     10GB
registry.hzdr.de/proteinprompt/proteinprompt   latest        cd99d22f3cb2   18 months ago   1.93GB
nvidia/cuda                                    11.0.3-base   4f3e95122fab   23 months ago   122MB

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker run  -v $(pwd)/data:/data proteinprompt search -f data/input.fa -o data/output.csv
Unable to find image 'proteinprompt:latest' locally
docker: Error response from daemon: pull access denied for proteinprompt, repository does not exist or may require 'docker login': denied: requested access to the resource is denied.
See 'docker run --help'.

Do you know what is wrong?

Final approach,
docker pull registry.hzdr.de/proteinprompt/proteinprompt:latest

On what operating system? Please [edit] your question and add your OS, and what you
have tried so far (and how it failed) so we don't tell you things you already know don't work.

@M__ thanks, but i have new issue... At least I'll learn the docker

@M__ I pull it docker pull registry.hzdr.de/proteinprompt/proteinprompt:latest :)

Okay, I see ... right. It's an alpine docker container. It will not have apt-get, there is a different command for apt-get in alpine its not apt-get but I don't know what it is. Sorry thats my fault. I was thinking it would be using Ubuntu or Debian - but no its alpine. Alpine is smaller thats why it's being used - but has very little in it. It would have been alot easier to use Ubuntu/Debian, the authors should have included git in the container.

Anwer is RUN apk update thanks for support!

The docker file is here https://codebase.helmholtz.cloud/proteinprompt/ProteinPrompt/-/blob/main/Dockerfile
You need to download and install a copy of Docker (below).
These are the instructions how to build the container locally from the remote Docker file. Copy and paste this code into a file called Dockerfile. Here I placed it in a tmp directory, then simply (see note)
 docker build -t tmp/dockerfile .

Just to note the authors may have specific instructions for the build command.
How to install Docker Before this command is run you need to install docker here https://www.docker.com/products/docker-desktop/
After installation you need a special command to run Docker in sudo-less mode is:
sudo usermod -aG docker

The first command I gave will then build the Dockerfile. There quite a lot of commands BTW to run/execute a Dockerfile, but you can also ask more questions.

Why am I recommending this way?
Simple, because normally downloading docker is done via https://docker.com and you would issue a pull request. For some weird reason the authors haven't done this and I don't know what the outcome of an external download is. If it is built via the Dockerfile then Docker will recognise the Docker container. Once the container is built locally you can upload it to your own Docker site and delete it locally and access it any time via a pull request or a run request (will automate a pull request).
Docker is very cool BTW, worth learning.

Comments, following bug report
I don't think you've got git installed (apt-get install -y git). Once that is installed I would would try again. BUT FIRST I would delete the current build docker rmi dockerfile. If docker rmi ... fails (I don't think it will) then remove all containers before deleting the image with docker rm $(docker ps -qa). The use rmi and then install git apt-get install -y git and do the docker build command.
If this doesn't work (same error) ..... you need to alter the Dockerfile and build it again from the new file .. Its a bug by the authors simply add the line RUN apt-get update && apt-get install -y git ... This is because there's a chance that git needs to be on the inside of the container not the outside.
If all the above fails, I would download the entire image/container from the website the authors have loaded it. I don't know how to pull from a non-docker website, I've never done it. Note I'll place these comments in an updated post. I am honestly skeptical about the authors ability to write/develop Docker.


text:
2D diagram of intermolecular interactions in R
I've docked a large inorganic compound (receptor) with an organic molecule that isn't a protein (ligand). I wanted to see the interaction between them in a 2D diagram, but BIOVIA Discovery Studio (DS) doesn't work well with pdb files converted from pdbqt. And since the receptor isn't a protein or a fragment of DNA/RNA, DS also won't show the diagrams (the same with PLIP and Proteins Plus). I've also tried LigPlot, LeView and PoseView, but none of them worked as well. Does anyone know another  software that draws 2D diagrams, and, if not, is there any way to do this in R???

Hi @ginn are you sure R is the right platform for this job? Generally, this is where Python excels.

Oh... I hadn't thought of that yet  but thank you so much for your input. Would you know, perhaps, of any package that does this job...?

Beyond Bio.pdb https://biopython.org/docs/1.75/api/Bio.PDB.html and of course pymol - I dunno - I'm not a specialist in this area. There are members here who are experts, so maybe they'll respond.

I've tried PyMOL, but apparently you have to install the pose view plugin, which you need to be licensed to use...  But thank you so much once again!!! I'll try bio.pdb in the meantime ☺


text:
Is a human genome file with reference as NCBI37 the same as hg19?
This question was also asked on Biostars
In other words, are NCBI37 and hg19 synonymous? I can't find much information for what NCBI37 is in a human (not mouse) genome.
The VCF is from a lab sequencing human genomes, and they're no longer able to be reached out to. I can't provide a link to NCBI37 because I don't know what it is and I can't find much information on it.
Relevant line from VCF:
"##reference=NCBI37"
Here is part of the VCF header:
##fileformat=VCFv4.2
##source=apt-format-result:2.10.2
##reference=NCBI37
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1>
##contig=<ID=10>
##contig=<ID=11>
##contig=<ID=12>
##contig=<ID=19> [omitted some chromosomes for brevity]
##contig=<ID=2>
##contig=<ID=20>
##contig=<ID=MT>
##contig=<ID=X>
##contig=<ID=Y>

Most likely yes, they are the same.
However, "NCBI37" is an odd/unusual notation, because it doesn't indicate the species (human / mouse). Assuming it's a human genome (based on you comparing it to hg19 rather than mm19, or something else), then the correct name for the NCBI reference should be GRCh37.
As it's a different name, we can only make guesses as to what the people who generated the VCF meant. You'll need to talk directly to them in order to be sure.


text:
How important are the homozygous variants that get unnecessarily deleted using liftover?
I'm referring to the text described here:

These tools [NCBI remap, CrossMap] operate only on the sites present in an input VCF, and return the representation of those sites in a new genome assembly. This does not capture all variation, however. Consider an individual who has sequence reads that indicate they match the GRCh37 reference genome assembly at position GRCh37.chr1:169,519,049 (i.e. the individual's genotype is T/T). Because the individual is homozygous reference at that site, there will be no variation present in their VCF file created on GRCh37. However, the analogous position on the updated GRCh38 reference genome assembly, position GRCh38.chr1:169,549,811, has the reference base C. Consequently, if the individual's read data were analyzed on GRCh38, they would be identified as homozygous for a C->T SNP. Because this site is not present in the input GRCh37 VCF, it is never added when creating a GRCh38 VCF by these other tools.

How important is this? Do you think it is relevant to any types of analysis?

I'm mostly just wondering if the bioinformatics community thinks the added benefit is large enough to be relevant. I haven't encountered a need for this in a project.

All this is saying is that the tool doesn't lift over variants which aren't present in the vcf. I can't imagine this being a big deal in any analysis I've done and it's exactly what I'd expect the tool to do. There's no rule that invariant sites (with respect to the reference) can't be in a vcf, so if you are concerned then just retain the variants in the vcf file at the calling stage


text:
How many protein sequences do we usually have for each species?
My question is exactly the one in the title. I add a particular example for more clarity. If I want to study the sequence of protein IkBα, do I find only one sequence or is there a database that contains many sequences from many human host species ?
EDIT: The analysis that I am going to do is to use NLP models to apply it to protein sequences to extrapolate, for instance, protein function knowledge.
Here is a reference of a research paper Learning the protein language: Evolution, structure, and function that considers this application by taking into account the problem of low number of sequences data.
So I think I should have at least 10 protein sequences.
Thank you in advance.

@M__ of course, I am going to read Maximillian's response and to answer !

Just to add to Maximilian Press's answer. I happen to know about human immunology and I think Maximilian Press has a good idea.
Genetic studies and immunological studies per se in cellular immunology  tend to be restricted within humans. There are limited exceptions to this namely chimpanzees and humanised mammals, most notably humanised mice. In cellular immunology these mammals retain immunologically distinct behaviour but are sufficiently similar to be useful models.
The protein in question is a B-cell repressor, thats more complicated ...
However, using trans-species homology of the protein of interest for

chimpanzees
humanised mice
regular mice

.. represents a good diversity of proteins to use what appears to be a deep-learning function prediction if this was a traditional cellular immunological study. Any immunologist looking at the biological rationale would recognise the model systems deployed.
Method: Blast Generating diversity for given protein is easy, except for humans. Just use online Blast and restrict the search to the species of interest and 20000 hits (?I forget the upper limit) is the maximum and select for more diverse genes. This would work for chimps and (possibly) mice. Thus blast a human gene and restrict the output to the species of interest. Human gene diversity - this much bigger than an online blast could be targeted in accordance with human genetic diversity. The cellular immune system for humans does not conform to background diversity in general. Some geographic human populations do have a HLA predominance but there are no absolute association.
The final option is local blasting the human database for that gene, that is doable and objectively used to derive representative diversity - for this purpose is complicated and would advise against it.

thanks @M__ . Even if you define your answer as an "addition" to the first one, it was what I was looking for (a general explanation about the database).

Thanks for your feedback @HelpNeederStudent

I'm going to start by assuming that you don't actually want 10 protein sequences from the same species, even though you write this in your comment on M__'s (deleted) answer. This will not be a very informative sample of variation in a given protein, except for extreme examples like antibodies.
To get a good sample of variation you almost certainly want to look across species. The paper that you link to does this by using the UniRef database, they then cluster proteins at the 90% similarity threshold. That is well beyond the variation that you'd expect within a species. See their section "Training datasets":

We train our masked language models on a large corpus of protein sequences, UniRef90 (Suzek et al., 2007), retrieved in July 2018. This dataset contains 76,215,872 protein sequences filtered to 90% sequence identity.

You could use an approach like that in the paper to get similar-ish clusters of proteins whose variation you can then analyze. You can use this, but I am somewhat suspicious of that, without reading the paper in too much depth, because we know rather a lot about orthology/paralogy and how it affects function in these cases. Especially in these very elaborate not very introspectable models I worry about what they will do with inadequately curated data.
So instead I'd probably want to first identify the orthology groups (gene families) of interest. These are generally manually curated groupings of proteins.
Once you have the orthology groups that you care about, you can use a database such as KEGG or OMA to identify the specific orthologs (genes) that belong to the group. Then you will have your dataset.

Thank you Maximilian for your explanation about the paper. However, if one wants to search for motifs, the dataset should be 10 protein sequences from the same species, right ? Even though, as @M___ has explained, "the cellular immune system for humans does not conform to background diversity in general. Some geographic human populations do have a HLA predominance but there are no absolute association".

@HelpNeederStudent I am not sure why you would need them to be from the same species. You will have to clarify that. The whole idea of a gene family or orthology group is that it's a bunch of curated genes/proteins that have the same essential biological function. The people in your paper certainly aren't concerned with species of origin! You could still theoretically do it within a species, I'm just skeptical that you'll find a level of variation that you can do any work with (my professional though non-expert opinion).

Thank you very much !


text:
Singularity and Illumina's Nirvana
I am trying to using Illumina's Nirvana however I can't seem to get it to work despite having the right Singularity path.

I pulled Nirvana down using

singularity pull docker://annotation/nirvana:3.14

I found the path to the .dll using the following command:

singularity exec nirvana_3.14.sif /bin/bash -c "cd /opt/nirvana/ && ls"

Tried to run the downloader using:

## Define path to singularity image:
sif_container="/path/to/nirvana_3.14.sif"

## Path in the singularity file:
sif_path="/opt/nirvana/Downloader.dll"

## Data to be placed in
Data="path/to/Nirvana/GRCh38/"

#_____________________________________________________________________________________________________________________________________

singularity exec "$sif_container" dotnet "$sif_path" --ga GRCh38 --out "$Data"

yet it keeps saying ERROR: The top-level output directory directory was not specified. Please use the --out parameter.. I have checked the paths and they are correct and I have specified the data path yet it does not seem to recognize it? I have even moved around the container and the path variables to no avail?

From your previous questions, it is apparent that your HPC has an odd singularity setup. You could ask you sysadmin about how to navigate it.
EDIT
After working with OP on chat (and based on previous experience with them), it looks like their sysadmin has mount hostfs enabled in the singularity.conf file, which messes with every container every time. I guess it's easier than letting people bind their directories of choice but anyway, the way to circumvent that is to use --no-mount hostfs
@Indira, always run singularity [exec|shell|run] --no-mount hostfs instead of plain singularity [exec|shell|run]. All containers will then run well and if they error out, it's probably because either the container internally encountered an error or it was unable to access an external resource. What you're running into now is a mix that obscures the real problems and creates a no-win scenario for you.


text:
Make a new column in R or bash
I have antibody data the first column is heavy chain sequence and the second column is the light chain
I want to create a third column by combining two columns like
If CARGRGIAVAGTTGFDYW is a heavy chain and CQQYNSYSLTF a light chain, in third column I want CARGRGIAVAGTTGFDYW_CQQYNSYSLTF
if a heavy or light chain is absent, I want something like NA_CQQYNSYSLTF or CARGRGIAVAGTTGFDYW_NA depends of absent part
CDR3.light.aa   CDR3.heavy.aa
-              CAGDLSVGGVLAGGVPHLRHFDPW
-              CPETPRYCSGGFCYGEFDNW

CSSYAGSNTFGLF   -
CSSYAGSNTGLF    -

So far I have done this
a$a <- paste(a$CDR3.heavy.aa, a$CDR3.light.aa, sep= "_")

But I noticed a lot of mess and manual works
Any help please?

What did you try so far?

Thank you, I have added my code in the post

Your approach was correct, you just have to transform your dataframe first to replace the character "-" by NA :
library(stringr)
library(dplyr)
df=data.frame(CDR3.light.aa=c("-","-","CSSYAGSNTFGLF",'CSSYAGSNTGLF'),
              CDR3.heavy.aa=c("CAGDLSVGGVLAGGVPHLRHFDPW","-","CPETPRYCSGGFCYGEFDNW","-"))

df=df %>%
  mutate_all(funs(str_replace(., "-", "NA")))

df$a <- paste(df$CDR3.heavy.aa, df$CDR3.light.aa, sep= "_")

 df
  CDR3.light.aa            CDR3.heavy.aa                                  a
1            NA CAGDLSVGGVLAGGVPHLRHFDPW        CAGDLSVGGVLAGGVPHLRHFDPW_NA
2            NA                       NA                              NA_NA
3 CSSYAGSNTFGLF     CPETPRYCSGGFCYGEFDNW CPETPRYCSGGFCYGEFDNW_CSSYAGSNTFGLF
4  CSSYAGSNTGLF                       NA                    NA_CSSYAGSNTGLF


text:
Singularity And MacOS
I am trying to install Singularity on my MacOS but it doesn't seem to pick up which singularity or singularity --version:
Main Attempt:

Install using Homebrew (and install Home-brew):

/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"
brew install --cask --no-quarantine singularity
singularity --version

resulting in
==> Downloading https://bitbucket.org/router_gray/singularityviewer/downloads/Singularity_Alpha_1_8_9_8709_x86_64.dmg
Already downloaded: /Users/Username/Library/Caches/Homebrew/downloads/d3015dad9cc2d3ef25603a68114f3caf63cb0dd90e03f159405999262392baec--Singularity_Alpha_1_8_9_8709_x86_64.dmg
==> Uninstalling Cask singularity
==> Backing App 'SingularityAlpha.app' up to '/opt/homebrew/Caskroom/singularity/1.8.9.8709/SingularityAlpha.app'
==> Removing App '/Applications/SingularityAlpha.app'
==> Purging files for version 1.8.9.8709 of Cask singularity
==> Installing Cask singularity
==> Moving App 'SingularityAlpha.app' to '/Applications/SingularityAlpha.app'

Then test using singularity --version with the result
singularity not found

Attempt 2:

Maybe it wants singularity spelled with a capital Singularity instead of lowercase:

which Singularity

resulting in:
Singularity not found

Attempt 3:

Export the path to the bash_profile and the source it

export PATH="/Applications/SingularityAlpha.app/Contents/MacOS:$PATH"

However, this resulted in the same issue when singularity --version was run. It could not recognize that the singularity was installed.
Attempt 4:

Explore inside the Singularity app which is installed under '/Applications/SingularityAlpha.app'

The path /Applications/SingularityAlpha.app/Contents/MacOS lead to a file called SingularityAlpha.

Exporting this path in the bash profile and sourcing still lead to the same issue when trying to see if it even recoginzed the singularity using singularity --version.

I have also tried reinstalling it and it still has the same issue. Is there any suggestion for how to install/use Singularity locally as many packages will have to be tested locally on my machine before running it on my HPC?

Answer: Currently there is not a good way for Singularity to be installed on MacOS. Perhaps, there will be in the future.

I've tried and failed to install Singularity on macOS - it needs VirtualBox  or some such virtualization software and a command line utility that makes virtual containers - ultimately there were so many problems that I just gave up. Docker is far easier to install on macOS and there is not a lot of difference IMO. You can use docker to build and test images and then pull them to a place where you have access to singularity and test singularity specific stuff there.
I think Singularity is a lot easier to install on Linux. See here: https://docs.sylabs.io/guides/3.0/user-guide/installation.html#mac

EDIT-2
Singularity-Alpha is a Second Life video game related cask (https://github.com/singularity-viewer/SingularityViewer), it is not the Singularity software. Don't waste your time on it, I just tried to test it and realized I'd already tested and discarded it as irrelevant.

EDIT: To address your homebrew issue, note that you're installing a cask and it's moving an App to /Applications/. Just typing in commands that should work is not the same as proper testing. Try opening the SingularityAlpha application to see what that does. I'm not saying it'll work, just that you can test it better. Your attempts only show a lack of understanding in how macOS Applications work and how you can(not) access them from the command line or add locations to $PATH like you can on Linux.


text:
overrepresentation test, between transcriptome and candidates sequences obtained from the transcriptome
For an analysis of my data, I have a transcriptome and a list of sequences obtained from the transcriptome. I would like to perform a functional enrichment analysis. I have annotated both sets of data using eggnog mapper. Currently, I want to perform a test between the two functional annotations, specifically COGs (Clusters of Orthologous Groups). I have tried using the R code https://yulab-smu.top/biomedical-knowledge-mining-book/enrichment-overview.html#gsea-algorithm
with clusterProfiler, but it seems that it may not work. With which tools or code can I perform this test, please?
exemple somme of my data :


text:
Singularity, AnnotSV and HPO
I am trying to use AnnotSV to annotate some Manta outputs in the following code with the HPO term changed:
Where:

input_dir is the input directory
singularity_container is the path to the AnnotSV singularity container

# Loop through the files in the input directory
for file in "$input_dir"/*.vcf; do
    # Get the base name of the file without extension
    filename=$(basename "$file" .vcf)

    # Create the output subdirectory within the input directory
    output_dir="$input_dir/AnnotSV_Batch_Run"
mkdir -p "$output_dir"

    # Set the output file path with "_AnnotSV" appended
    output_file="$output_dir/${filename}_AnnotSV_"

    # Run AnnotSV via Singularity with the base name as output prefix
    singularity exec "$singularity_container" AnnotSV -SVinputFile "$file" -genomeBuild GRCh38 --hpo "HP:0000001"  -outputFile "$output_file"
done

The HPO term I am using is in the jax.org site but when I run it - I get the most peculiar error  option not known. but the HPO term is specified in the AnnotSV manual on page 49 (according to the number on the bottom). As for trying I tried with --HPOterm and hpoterm I tried with and without the : and with and without quotes as well. I am not sure what the issue is but I cannot seem to find the solution. Has anyone perchance used this software before with HPO terms as it seems not to be accepting the input I give it?

The manual says -hpo and you're using --hpo. Could that be the problem?

Now I am getting "/usr/local/share/AnnotSV/Annotations_Human/Genes/GRCh38/refGene.txt.gz" file doesn't exist
Please check your install - Exit with error.?

Try adding a --no-mount flag immediately after singularity exec in the command.

The --no-mount gives the error FATAL:   "--no-mount": executable file not found in $PATH - I also tried with usr/local/bin/AnnotSV and it did not work

Oops, my bad. I did not read the documentation properly, sorry. Try singularity exec --no-mount hostfs AnnotSV

This error is persistent and I am not sure how to get around it.

Is this the page from which you're getting the singularity image link: https://lbgi.fr/AnnotSV/downloads

Try downloading the annotations file, uncompressing it and providing it as an input to the -annotationsDir parameter. The base image, which is ~400 MB seems to contain only the executable and all the images I see on docker hub are at least 4 gigs, which makes me think they include some sort of resource file that the base image does not have.

Yes the singularity is from the link

I believe I am doing this incorrectly as I give it the folder of the phenotype_to_genes.txt that I got from HPO but it seems like it still gives command not found for a lot this is only part of the error ```/usr/local/etc/AnnotSV/configfile: line 9: -annotationMode:: command not found
/usr/local/etc/AnnotSV/configfile: line 11: -benignAF:: command not found
/usr/local/etc/AnnotSV/configfile: line 13: -candidateGenesFiltering:: command not found
/usr/local/etc/AnnotSV/configfile: line 15: -genomeBuild:: command not found

The folder I gave was under the annotationsDir flag but the error persists.

Something is interfering with the shell parsing your command. Add an echo in front of the singularity command, then add a set -x before the for loop begins. Run the script and paste the output on pastebin, then link to the pastebin content here.

What is pastebin?

Did you Google it?

I did - it seems that AnnotSV doesn't work well with Singularity so this could not be done.

It seems that currently AnnotSV does not support HPO terms usage an Singularity.


text:
Help with running ATAC-seq using Encode pipeline
I am very new to ATAC-seq pipeline from Encode and trying to run on a HPC but not sure the steps after reading their instruction. https://github.com/ENCODE-DCC/atac-seq-pipeline
Here is my json file:
{
    "atac.title" : “atac)”,
    "atac.description" : “Encode”,

    "atac.pipeline_type" : "atac",
    "atac.align_only" : false,
    "atac.true_rep_only" : false,

    "atac.genome_tsv" : "https://storage.googleapis.com/encode-pipeline-genome-data/genome_tsv/v4/hg38.tsv",

    "atac.paired_end" : true,

    "atac.fastqs_wt_rep1_R1" : [ "rep1_R1_L1.fastq.gz" ],
    "atac.fastqs_wt_rep1_R2" : [ "rep1_R2_L1.fastq.gz"],
    "atac.fastqs_wt_rep2_R1" : [ "rep2_R1_L1.fastq.gz" ],
    "atac.fastqs_wt_rep2_R2" : [ "rep2_R2_L1.fastq.gz" ],
    "atac.fastqs_di_rep1_R1" : [ "rep1_R1_L1.fastq.gz" ],
    "atac.fastqs_di_rep1_R2" : [ "rep1_R2_L1.fastq.gz"],
    "atac.fastqs_di_rep2_R1" : [ "rep2_R1_L1.fastq.gz" ],
    "atac.fastqs_di_rep2_R2" : [ "rep2_R2_L1.fastq.gz" ],

    "atac.auto_detect_adapter" : true,

    "atac.multimapping" : 4
}

I also need input file so what it should be?

Input files
IMPORTANT: Our pipeline considers a replicate (rep) as a biological replicate. You can still define technical replicates for each bio replicate. Tech replicates will be merged together to make a single FASTQ for each bio replicate.
IMPORTANT: Our pipeline supports up to 10 bio replicates.
IMPORTANT: Our pipeline has cross-validation analyses (IDR/overlap) comparing every pair of all replicates. Number of tasks for such analyses will be like nC2. This number will be 45 for 10 bio replicates. It's recommended to keep number of replicates <= 4.
Pipeline can start from any of the following data types (FASTQ, BAM, NODUP_BAM and TAG-ALIGN).

Starting from FASTQs

Technical replicates for each bio-rep will be MERGED in the very early stage of the pipeline. Each read end R1 and R2 have separate arrays atac.fastqs_repX_R1 and atac.fastqs_repX_R2. Do not define R2 array for single-ended replicates.

Example of 3 paired-ended biological replicates and 2 technical replicates for each bio rep. Two technical replicates BIOREPX_TECHREP1.R1.fq.gz and BIOREPX_TECHREP2.R1.fq.gz for each bio replicate will be merged.

       {
           "atac.paired_end" : true,
           "atac.fastqs_rep1_R1" : ["BIOREP1_TECHREP1.R1.fq.gz", "BIOREP1_TECHREP2.R1.fq.gz"],
           "atac.fastqs_rep1_R2" : ["BIOREP1_TECHREP1.R2.fq.gz", "BIOREP1_TECHREP2.R2.fq.gz"],
           "atac.fastqs_rep2_R1" : ["BIOREP2_TECHREP1.R1.fq.gz", "BIOREP2_TECHREP2.R1.fq.gz"],
           "atac.fastqs_rep2_R2" : ["BIOREP2_TECHREP1.R2.fq.gz", "BIOREP2_TECHREP2.R2.fq.gz"],
           "atac.fastqs_rep3_R1" : ["BIOREP3_TECHREP1.R1.fq.gz", "BIOREP3_TECHREP2.R1.fq.gz"],
           "atac.fastqs_rep3_R2" : ["BIOREP3_TECHREP1.R2.fq.gz", "BIOREP3_TECHREP2.R2.fq.gz"]
       }

How should I modify this command:
caper hpc submit atac.wdl -i "${INPUT_JSON}" --conda --leader-job-name ANY_GOOD_LEADER_JOB_NAME

I appreciate your help !

This is very impressive - cross-validation cool. The thing I would check is this is can be run via conda, docker and singularity (you are using --conda). You need to know which one is recommended for your HPC.

Thank you . I usually run singularity.


text:
Reading same BAM file twice with htslib
I would like to read iteratively the same BAM file twice using htslib.
My option is to use hts_open, sam_hdr_read, then sam_read1 as much as needed to read the file until EOF.
Is there a possibility to "rewind", and use sam_read1 again from the beginning of the file, without closing the file, destroying the object, and opening it again?

EDIT
My aim is group the reads by the barcode.  The first pass would allocate the specific sizes first (the number of reads by barcode), then store the reads.
I know I could sort the reads first (using samtools sort -t BX), but I would like to skip this step using pre-allocation.

@M__ Thanks, I edited the question.  I hope it helps!

The best solution is to reopen the file. For a coordinate sorted file you have the option to jump back to the start of the first chromosome with index but that is more complex, less flexible and won't be faster.


text:
How to calculate the variant allelic fraction (VAF) from VEP output?
This question was also asked on Biostars
How can I calculate the variant allelic fraction from the VEP file for multiple human samples?
I used the ensembl-vep and it is human data.

Please show a) your VEP output and b) describe what variant allelic fraction  means. Thank you.

If you just have a VEP output (i.e. VCF) without raw genotypes then you cannot calculate the allelic fraction.

thanks user438383, how can calculate it if i have the vcf file?

@Joman See: https://www.biostars.org/p/147443/


text:
Residue numbering in biopython
Is there an out-of-the-box function in the Bio.SeqIO module (in biopython) that can print a protein sequence (in command line) in multiple lines along with the numbering of the residues (either left and right or in-between the lines)?
In particular I have a .fasta file which I can read into a SeqRecord object like so:
from Bio import SeqIO
record = SeqIO.read("protein.fasta", "fasta")

I would like to then print the Seq class with the numbering, not just as a string as shown below:
print(record.seq)

I have already tried:
SeqIO.write(record, sys.stdout, 'gb')

But it is not working for me. The error was ValueError: missing molecule_type in annotations. I googled it and it appears to be related to newer biopython versions (see https://github.com/linsalrob/PhiSpy/issues/43).
This on-line tool seems to do what I want, but one needs to provide it with an alphabet. How do I do that with code? Does a SeqRecord object have an alphabet?

@M__ is right - there's a molecule_type error when you try to write and this needs to be fixed. Here's my end to end approach:
The input file:
cat xyz.fa
>NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAA
PPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKT
CPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRN
TFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGR
DRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALEL
KDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

My (bad) code:
from Bio import SeqIO

with open("xyz.fa") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        record.annotations = { "molecule_type" : "protein" }
SeqIO.write(record, "xyz.gb", "gb") # moved two indentations to the left

The output:
cat xyz.gb

LOCUS       NP_000537.3              393 aa                     UNK 01-JAN-1980
DEFINITION  NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens].
ACCESSION   NP_000537
VERSION     NP_000537.3
KEYWORDS    .
SOURCE      .
  ORGANISM  .
            .
FEATURES             Location/Qualifiers
ORIGIN
        1 meepqsdpsv epplsqetfs dlwkllpenn vlsplpsqam ddlmlspddi eqwftedpgp
       61 deaprmpeaa ppvapapaap tpaapapaps wplsssvpsq ktyqgsygfr lgflhsgtak
      121 svtctyspal nkmfcqlakt cpvqlwvdst pppgtrvram aiykqsqhmt evvrrcphhe
      181 rcsdsdglap pqhlirvegn lrveylddrn tfrhsvvvpy eppevgsdct tihynymcns
      241 scmggmnrrp iltiitleds sgnllgrnsf evrvcacpgr drrteeenlr kkgephhelp
      301 pgstkralpn ntssspqpkk kpldgeyftl qirgrerfem frelnealel kdaqagkepg
      361 gsrahsshlk skkgqstsrh kklmfktegp dsd
//

I used the code from https://github.com/biopython/biopython/issues/3156 to fix the ValueError: missing molecule_type in annotations that I faced when I'd skipped the record.annotations = {"molecule_type" : "protein" } line.
There is definitely a better way to write the python code instead of using a loop or at least writing to a static file from within a loop, but that's not central to the problem here so I'm choosing to ignore it. I welcome any improvements to the code.

I've moved this two indentations to the left (one indentation left is okay). The example will work because it's just one sequence. Any more sequences will trigger overwriting.

Given the criteria required the closest solution is:
from Bio import SeqIO

seq_obj = ... # Biopython seq_obj
with open("mysequence.gb", "w") as fout:
    SeqIO.write(seq_obj, fout, "gb")

If this fails I would substitute gb for fasta and try again. If that fails, i.e. no fasta, there's a problem with the sequence object. If this succeeds but the gb fails, just supply a reproducible example. It will work.
I've used gb to successfully write APIs against NCBI Genbank and trap the meta-data. Genbank format is really cool.
Oh wait, sorry the answer is  ...
seq_obj for records:
    print(seq_obj.format("gb"))

From the comments "molecule type"...
Which are you using
for record in SeqIO.parse(filename, "fasta"):

or
for record in SeqIO.read(filename, "fasta"):

Which ever is used I'd to the other method.
Just to comment on @RamRS solution, thats cool I didn't know.
for record in SeqIO.parse(filename, "fasta"):
    record.annotations = {"molecule_type" : "protein" }

What is happening is the instance is being manually set. However, Biopython should detected the data type automatically. My suspicion is that read or parse one comprises and automated data type setter.
In any case, RamRS's solution is useful and will work. In past coding I would have read the object via Genbank.
The old alphabet method is deprecated and shouldn't be attempted.

For future reference RamRS's code is described here

I fixed it using this: with open("xyz.fa") as handle: for record in SeqIO.parse(handle, "fasta"): record.annotations = {"molecule_type" : "protein" } SeqIO.write(record, "xyz.gb", "gb") (excuse the bad python)

Code from: https://github.com/biopython/biopython/issues/3156

@RamRS Maybe possibly give Leptoceratops an upvote? This question has been a long haul for everyone, but it's a good answer.

Sure thing, they have been very cooperative and clear. I appreciate their contribution


text:
Gviz Coverage Plots
This question has also been asked on Biostars
I have two bam files from single-cell RNA sequencing mapped to the reference genome using CellRanger, I can view them in IGV and I have a particular region where the pattern of reads mapped to the reference genome are different between the two bam files but when I try the following code to plot trackplots with coverage, the resulting plot is not consistent with what I see in IGV. I am not very familiar with this, aren't I supposed to normalize the data as well? please help me.
txdb <-TxDb.Hsapiens.UCSC.hg38.knownGene
txTr <- GeneRegionTrack(txdb, chromosome = "chr1",
                        start = 153561336,  end = 153561564)
dTrack <- DataTrack(range = bam.file, genome = "hg38", type = "l",
                     name = "Coverage", window = -1,
                     chromosome = "chr1")

dTrack0 <- DataTrack(range = bam.file.0, genome = "hg38", type = "l",
                     name = "Coverage", window = -1,
                     chromosome = "chr1", col="red")

plotTracks(c(dTrack,dTrack0, txTr), from = 153561357, to = 153561585)

Welcome @A4747 its useful to describe what you are dealing with, e.g. humans, a model system and any further project details. Essentially the answer will be project specific. Anyway provided an appropriate summary paper.

I would recommend reading the following study by Zhao et al (2021).
The traditional simplest measure is CPM otherwise known as TPM. This is counts/transcripts per million. FPKM (fragments per kilobase of transcript per million fragments mapped) is used for metagenomics. I don't know its application in your context ... I don't know your context.
The discussions for the pros- and cons- of FPKM and TPM/CPM are discussed therein.
The authors in the study favour normalisation - in their context - not in general. The authors use the coefficient of variation, and intraclass correlation coefficient. I would see this as a form of standardisation. It's worth keeping an open mind because its a pivotal step in the calculation but all too easily overlooked area of investigation.
From RamRS's comments that may have changed already. In general statistics standardisation is a pivotal step, usually followed by the removal of outliers.

Unfortunately, neither TPM nor CPM can be used in any context to compare two samples. FPKM is obsolete IMO

FPKM is used in metagenomics (a lot, lot more than one species present). Recovering a full length gene/transcript rarely happens.

Thank you, I did not know that.

@M__ Thanks for your response. I agree that in this case TPM normalization may be required when displaying them in a trackplot. I think those types of normalization may not be accurate when performing differential expression analysis through DESeq. However, my particular problem is with the fact that what I see in IGV ends up being different from what I see when I plot the bam file in R through the script written above.


text:
How do I calculate differential expression for RNA-seq values with the "limma" package and the "ebayes" function?
So for context, I have a set of TPM values for multiple genes for different samples, and I need to calculate the differential expression for RNA-seq values. I've been using this website as a guide: https://raw.githubusercontent.com/ucdavis-bioinformatics-training/2019_March_UCSF_mRNAseq_Workshop/master/differential_expression/DE_Analysis.Rmd. The overall aim of this calculation is so that I can use PANDA to map TM motifs/combine gene expression data to create networks representing interactions between transcription factors and genes, and GSEA to analyze genes ranked by fold change or differential expression p-value.
So what I need to do is match one male to one female of the same age group (which I already did prior; the number of males outnumbered the number of females so I filtered the number of males so that there would be an equal number, and of the same age distribution). But it seems that I need to derive two factors and create a new variable "group" that combines factor1 and factor2. In this case, it would be the gender I would use as the group factors, right?

My advice is open a separate question on the topic of calculating DEG cut-offs on this site. I personally cannot  answer the first question (I don't know the project, specifically what the DEG is measuring and how pairwise DEGs are avoided), but I can present the rationale of calculating DEG cut-offs (i.e. different approaches used and why). Again this is the second part of the question and the SE (Stackexchange) approach is one technical request per question. In this instance the answer  is bit complicated (if stats aren't you thing) and is not a one line response in comments.

Okay, but do you have a recommendation on the specific number I use as a cutoff?

Which is..? It wasn’t clear to me if you had a number in mind or you were telling me to just eyeball it.

And how do you derive this?

Okay, I did what you asked


text:
Mapping single cell data with annotation
Entry level Bioinformatician here.
I am trying to map my single cell data with standard annotation file or assign my Seurat cluster to cell types. I have used the singleR and celldex package to map it but my problem is the single cell only consist of endothelial cells. I haven't found paper that only consist of rats endothelial cells. The single cell data is for rats.  I am stuck. What is the best thing that could be done here?

Given that it sounds like you already have a Seurat object, can you please show a UMAP of the cells that is coloured based on clusters identified in FindClusters? This will help us to understand the heterogeneity of the cells you are looking at.

It seems like you can annotate all your cells as endothelial cells, problem solved ;) But seriously, not clear what your expected cell types are (subtypes of endothelial cells perhaps) or the tissue type. If you can't find a suitable reference dataset for the reference mapping approach, then you can always resort to the 'old fashioned' way of unsupervised clustering (e.g leiden) plus checking marker gene expression to annotate the clusters.
An alternative approach, could be finding a suitable human/mouse reference and  then converting the gene ids to rat (ideally via 1:1 ortholog mapping) but this could get messy, so I suggest the former simpler way, at least as a starting point

Thanks, this seems to clarify the question. Is this what you mean @appleball?

Hello. thanks for the comments and concerns.
For sequencing the data, only endothelial cells were selected so the expected cell types would be the subtypes of endothelial cells, venous capillary or arterial. I have not had any luck in finding the annotation file for the endothelial cells.

Note that many cell selection methods are imperfect, so for your first pass you should start with confirming their identity by looking at the expression of markers of endothelial cells, and markers of other cells types (pericytes, muscle, ...).


text:
Phylogenetic tree from a large FASTA file
I want to create a phylogenetic tree using a 3.8MB FASTA file with ambiguous code. How can I perform multiple sequence alignment and the construction of phylogenetic tree with these ambiguous codes?
Could you please suggest some programs or R packages to complete this task?
How does MEGA11 software handle ambiguous codes? Will the alignment consider them to be mismatches?

Welcome @user17615. Firstly, I would recommend starting with a pilot study and then scale up. What is the genetic distances or taxonomic groups between the sequences, e.g. species, genera, families, orders (thats quite important)? Is it essential to do this in R, or what the requirements for R, e.g. is that the language you work in? Finally, how is your bash/Linux command line?

MegaX is poor for this purpose, it handles sequence ambiguities. The file sizes for complete phylogenetic analysis within R is going to be very tough, but its doable using "soft statistical approaches" for phylogenetic methods and can be processed via the ape library.