semantic similarity measurement for cell line ontologies

About the metric you use, there are several reviews on the topic, if I recall well I think that the information content based are the best ones, of those, I only recognized the JiangConrathSemSim listed on the webpage.

However, if you want to use a semantic similarity measure I would use the BMA method, which is more comprehensive than other methods. But the BMA method is only a way to aggregate several metrics of ontologies into a single numeric value. I am not confident that FastSemSim implements correctly those methods.

Note, that according to the license, you should contact the authors before using it:

  You must contact us before using it or its results or any application/tool based on it in any published work.

I am not sure that the cell lines fit in a ontolgogy. Because ontologies are defined as directed acyclic graph (DAG) and there are some cycles that happens in cells like the induced pluripotent stem cells show. However they might be right if you just compare the cell lines and not the state of the cells on those lines (use it only if you want to compare the lineage of your cells lines).


Building STAR Genome Index for nanopore RNA sequencing

The parameter is used to determine how much sequence STAR indexes on each side of a splice junction to improve its alignment accuracy.  For very long reads, this may not be ideal. I am not sure if STAR is capable of including multiple splice junctions since a long read is more than likely to span more than one.

It may be worthwhile to consider aligning your reads directly to transcripts and taking the reads that align poorly and using them to look for novel transcripts.


Trimmomatic trims most of the reverse strand

Maybe you need to remove the unpaired reads:
trimmomatic PE -phred33 -threads 2 \
    read1.fastq.gz read2.fastq.gz \
    read1_trimmed.fastq.gz read1_unique.fastq.gz \
    read2_trimmed.fastq.gz read2_unique.fastq.gz \
    LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15

Then you can try with bwa.


Cutting the sequence into several sequences with the information of a dataframe

We wouldn't write the code for you from scratch (thats your job), but its an easy problem because in Python any string is an "array" or "dictionary". Thus

bacteria1 = 'GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDVIEEXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
print (bacteria1[0:120])

This will define the first gene in your list, remembering the array starts at 0 whilst your gene starts at position 1 (so you need to subtract 1). You simply parse the input list so the gene lengths form a dictionary and loop through the dictionary of gene lengths ensuring the last length is added to the current one and then issuing a write command.

You don't really need biopython in this instance, unless it will parse the gene length file you need to open (I'm not aware it can).

BTW the format you supplied isn't fasta it looks like a Biopython statement, again in this problem I wouldn't bother because Python is perfectly capable of splicing a string.


A PDB file has strict spacing rules as documented here.
The coordinates should be:
31 - 38       Real(8.3)      x             Orthogonal coordinates for X.
39 - 46       Real(8.3)      y             Orthogonal coordinates for Y.
47 - 54       Real(8.3)      z             Orthogonal coordinates for Z.

So trying this on your line:
'HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000'[30:39]

I get 9.051 0., which is not a number.
The whole line is formatted incorrectly however.
So as you say you are using Autodock why not pass a mol2 to prepare_ligand4.py instead which is an acceptable format?
Alternatively, you could convert your multi-entry sdf to pdb files via RDKit in Python:
suppl = Chem.SDMolSupplier(sdfile)
for mol in suppl:
      # mol = Chem.AddHs(mol) # has hydrogens?
      # AllChem.EmbedMolecule(mol) # is 3D?
      AllChem.ComputeGasteigerCharges(m) # sdfs do not have partial charges
      Chem.MolToPDBFile(m, f"{mol.GetProp('_Name')}.pdb")


Handling sample identity in aggregated 10x libraries?

You have to keep track of the file order you used for aggr.  The suffix number represents which represents the aggregated sample. You could store this information in a text file and load it into R independently of the Read10X function and combine it with the cell names to get their sample names. From the page you link:

  Gem Groups

  This is an integer that is appended to each barcode in the
  gene-barcode matrix. For example, AGACCATTGAGACTTA-1 and
  AGACCATTGAGACTTA-2 are distinct cell barcodes from different
  libraries, despite having the same nucleotide sequence.

  The numbering of the GEM groups will reflect the order that the
  libraries were provided in the Aggregation CSV.


MSA (protein) with biopython or something else?

In my personal experience, MUSCLE is the easiest program to use in conjunction with Biopython. Biopython features a command line wrapper for this program, which makes it very easy to use. Make sure to download the appropriate MUSCLE program from drive5 and save it somewhere. E.g., if you are using Jupyter in Linux:
!wget https://www.drive5.com/muscle/downloads3.8.31/muscle3.8.31_i86linux64.tar.gz
!tar -xzvf muscle3.8.31_i86linux64.tar.gz
!cp muscle3.8.31_i86linux64 /usr/local/bin
!chmod 755 /usr/local/bin/muscle3.8.31_i86linux64

Then, you can run MUSCLE like so:
from Bio.Align.Applications import MuscleCommandline
def runMUSCLE(infile, outfile):
    muscle_exe = r"/usr/local/bin/muscle3.8.31_i86linux64" #Here is where we installed MUSCLE
    muscle_cline = MuscleCommandline(muscle_exe,
                                     input=infile,
                                     out=outfile,
                                     clwstrict=True #Output in clustal format (more visually pleasing), otherwise the output is in FASTA. Whichever you need.
                                    )
    muscle_cline()

You don't need to specify that your sequence is amino acidic, however keep in mind that the input file must be in FASTA format.


Are there computational tools to extract features of DNA sequences?

If you want to do this with biopython, the SeqUtils package could be a solution.


Pulling out a certain gene in a volcano plot

As you did for labeling genes with an adjusted p value below 0.05, you can subset your dataset for keeping only rows corresponding to "Casp14":

library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),
                  aes(label = Feature.Name), color = "red")

With this code, you should see now the labeling of the gene of interest (Casp14) in red on your volcano plot.


The easiest solution is finditer
use re

# Compile you search term where [A|T|G|C] means A or G or C or T
p = re.compile("[A|T|G|C]GG")

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())

** Example Output**
(0, 2)
(22, 24)
(29, 31)

Easy does it. Its worth looking at the documentation for re here, https://docs.python.org/3/howto/regex.html for alternative approaches using match or search below. finditer is definitely the solution in context.

I support your coding approach because it provides lots of flexibility. For example, you might want to target GG or TTT around a given locus and you can easily modify the code to accomodate this. Someone elses code would be less flexible.
In summary, I kinda like your substrings and do like the count method. I'm not a fan of your == in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via re.
The basic issues are,

Biopython would be used via AlignIO or SeqIO. This will allow any given format to be analysed. Yours approach is okay but prone to bugs (below). [probably a minor issue]
Your sequence is in some sort of weird format, like extended phylip format and it is not fasta (that wouldn't work). It needs verifying because removing the first 65 characters might remove bits of sequence [probably a minor issue]
Chopping of -2 from the length ... I dunno I'd check out -1 [minor]
C-style matching. I don't like this. I would definitely use re. One solution uses match or search, but the solution above uses finditer, Reg-ex will allow to search for NGG and TTTN, which you are not doing
I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style == would match this because thats 2 and 3 bp respectively. Just thoughts. You will have 'edge effects' at the beginning and end of sequences, e.g. TTT 3' will match but that isn't TTTN.
I like your count method, but I think finditer is better because it gives you location information
The dictionary is sort of fine. Personally, I would dump the dictionary into a pandas dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.
the if and elif would be placed in a subroutine, this is good practice to convert to an object orientated approach.


How can I subset WGS data to the level of WES variants?

bcftools would be my choice, I'm sure bedtools could do the trick, too.
Something along this
bcftools view --regions-file

or  --targets, --regions-file might require a tabix index.


It looks like you were using an old annotation. The problematic IDs you posted existed in the GRCh37 annotations, but don't in the most recent GRCh38 annotation. For that reason they were excluded. The IDs that have - as symbols don't have associated symbols, but are present in the database.

To use an archived version in biomart:

mart = useDataset("hsapiens_gene_ensembl", useEnsembl(biomart="ensembl", version=84))

That's an example for release 84.


bcftools filtering all files in a directory

I assume your shell is looking for the variable $f_out, but can't find it because you did not define it. It probably generated a hidden file with the name .vcf in your output directory.

Instead, use the following: -o ${f}_out.vcf. This will "protect/limit" the bash variable to just $f rather than $f_out.

When encountering issues like this it's worth putting echo in your loop right before your command, in this case bcftools (adding quotes whenever necessary). That would show you the command as how it is going to be executed and would show you, in this case, that the output name is incorrect.


clusterProfiler, groupGO: meaning of GeneRatio

As you seem to be learning, the clusterProfiler package is a bit "special". It doesn't filter the results you're seeing for relevance, it returns everything. The GeneRatio is the number of genes in your input list associated with the given GO term (e.g., 5) and then a / and the total number of input genes. This is generally not that useful. Since many things aren't documented, you may find the source code helpful.


Modeling number of reads mapped to a gene

This is quite a similar question to this one, but I notice that you're specifically asking about a probability distribution rather than differential expression. That sounds like gene dispersion and/or shot noise; a good reference for me on that has been the documentation for DESeq2:
http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The gene dispersion issue is discussed in more detail in their 2014 paper:
https://doi.org/10.1186/s13059-014-0550-8

We model read counts K ij as following a negative binomial distribution (sometimes also called a gamma-Poisson distribution) with mean μ ij and dispersion α i . The mean is taken as a quantity q ij, proportional to the concentration of cDNA fragments from the gene in the sample, scaled by a normalization factor s ij , i.e., μ ij =s ij q ij . For many applications, the same constant s j can be used for all genes in a sample, which then accounts for differences in sequencing depth between samples. To estimate these size factors, the DESeq2 package offers the median-of-ratios method already used in DESeq. However, it can be advantageous to calculate gene-specific normalization factors s ij to account for further sources of technical biases such as differing dependence on GC content, gene length or the like, using published methods, and these can be supplied instead.


I think you are looking for gsub function:

a$Tumor_Sample_Barcode <- gsub("A","res",a$Tumor_Sample_Barcode)

If you prefer, you can use the str_replace_all function from stringr package. But it tooks a little bit more time than gsub:

Dummy data

seq = as.factor(paste0("A",1:21))
library(data.table)
dt <- data.table(Barcode = rep(seq,10000),
                 Value = rnorm(210000))

Test of gsub vs str_replace_all

> system.time(dt$Barcode <- gsub("A","res",dt$Barcode))
   user  system elapsed
  0.077   0.000   0.077

> system.time(dt$Barcode <- str_replace_all(dt$Barcode,"A","res"))
   user  system elapsed
  0.123   0.004   0.126
```


Trouble using biomaRt to retrieve hgnc symbols from Ensembl transcript ids

You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862":

You can do this with one more line:

g <- gsub("\\..*", "", rownames(txi.kallisto$counts))
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart = mart))
##  hgnc_symbol chromosome_name ensembl_transcript_id
##1          NA              14    ENSMUST00000177564
##2          NA               6    ENSMUST00000178537
##3          NA               6    ENSMUST00000178862
##4          NA              12    ENSMUST00000179520
##5          NA              14    ENSMUST00000179664
##6          NA              14    ENSMUST00000196221

With these transcripts it doesn't seem that there is an equivalent hgnc_symbol. But is recognized, otherwise they are not recognized:

(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = rownames(txi.kallisto$counts), mart = mart))
## [1] hgnc_symbol           chromosome_name       ensembl_transcript_id
##<0 rows> (or 0-length row.names)

However you might be interested in doing something like this:

human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = "ensembl_transcript_id",
       filters = "ensembl_transcript_id", values = g,mart = mouse,
      attributesL = "hgnc_symbol", martL = human)
## The query to the BioMart webservice returned an invalid result: the number of
## columns in the result table does not equal the number of
## attributes in the query. Please report this to the mailing list.


How do I create a for loop to filter through different FDR values?

For loops in most languages can be done in a variety of ways. What you have currently coded, loops over values 1 to 3 by using

for (i in 1:length(fdr.list)){
...
}

However, you probably want to use a different form and loop over each value in that vector using something more like:

for (fdr in fdr.nums.vector){
...
}

This way you have convenient access to the FDR value in each loop via the fdr variable. Which you can then use in constructing your output filenames. Thus your code would become something like:

fdr.nums.vector = c(5,10,20)

for (fdr in fdr.nums.vector){
  atac.glmtop.males.opening.fdr[fdr] <- topTags(
    atac.glmfit.males,
    n = nrow(atac.dge$counts),
    sort.by = "none",
    adjust.method = "fdr"
  )$table %>%
  tibble::rownames_to_column("peakco") %>%
  filter(logFC < 0) %>%
  filter(FDR < (fdr/100)) %>%
  merge(.,annotated.atac.pbmc.narrow, by="peakco" ) %>%
  tibble::column_to_rownames("peakco")

  atac.glmtop.males.opening.fdr[fdr]

  write.table(
    atac.glmtop.males.fdr[fdr].opening[,c("chr", "start","end")],
    file = paste0("~/Downloads/sexM.age.groupHO.fdr",fdr, ".opening.bed"),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
  write.table(
    atac.glmtop.males.fdr[fdr].opening["GeneName"],
    file = "~/Downloads/sexM.age.groupHO.fdr"+ fdr + ".opening.GeneNames.txt",
    sep = "\t",
    quote = FALSE,
    row.names = FALSE,
    col.names = FALSE
  )
}

Advancing From for Loops to apply

As a beginner, for loops seem quite logical and easy to understand. However, R has a more powerful way to achieve the same thing through what is called "vectorisation" using the apply function or one of it's relatives lapply, sapply etc. Once you have understood and mastered for loops, you should take a look at apply, especially if you ever need to do a for loop many times or you find yourself doing nested for loops. R is notoriously slow with for loops and apply will likely be many times faster (potentially orders of magnitude). Here are some links for when you are ready to advance from for loops to vectorisation.

https://www.r-bloggers.com/how-to-write-the-first-for-loop-in-r/
https://www.r-bloggers.com/using-apply-sapply-lapply-in-r/


convert tRNA coordinates Result to Fasta format

You can select “Output BED format” on the tRNAscan-SE server:

Given the output BED format and your input reference in FASTA format, you can easily generate a FASTA file of tRNA sequences via bedtools getfasta:

bedtools getfasta -name -s -fo ‹output.fasta› -fi ‹reference.fasta› -bed ‹tRNAscan-SE-output.bed›


significant SNPs to annotated candidate genes

You can use BEDOPS' closest-features tool. You only need to have the SNPs coordinates in .bed format and a .bed file with gene locations. BEDOPS will look for the closest genes upstream and downstream of your SNPs. The command will look something like this:

closest-features yourSNPs.bed yourAnnotation.bed

Then, you can obtain your functional annotations through different tools such as biomart.


how to merge more than two sample in Seurat?

You have fed arguments to the MergeSeurat() function that it does not expect. In terms of objects, MergeSeurat() accepts only 2 arguments object1 and object2, please do ?MergeSeurat and see yourself.

Seurat is probably the best documented single cell package, your exact question has its own dedicated vignette: https://satijalab.org/seurat/v3.0/merge_vignette.html

As the vignette above is for Seurat v3 and in case you would like to use Seurat v2:
You can create a named list of Seurat objects to be merged and then use reduce() of the purrr package.

# create a named list of seurat objects to be merged
# code not shown

# optional but probably a good idea
# rename cells using object names as prefix
for (i in names(seurat_object_list)) {
  seurat_object_list[[i]] <- RenameCells(seurat_object_list[[i]],
                                         add.cell.id = i)
}

# merge all the objects in the list
merged_combined <- reduce(seurat_object_list,
                          MergeSeurat,
                          do.normalize = FALSE)


The linear gap penalty has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base.
The simplest case is assigning each unit (= each missing base) a penalty of 1:

  g(k) = k

Or more generally, if you want to increase the penalty per base to 2, 3 or n:

  g(k) = n * k
  with k = length of the gap, and n = cost of 1 missing base

The affine gap penalty has a cost that follows -as its name conveniently states- an affine function.
The idea behind it, is that any gap in an alignment has a cost with two components:

the "gap opening" penalty: it is a cost payed only once per gap, basically just for existing.
the "gap extension" penalty: equivalent to the linear gap penalty, it costs 1 unit for each consecutive missing base.

The general formula for the cost is then:

  g(k) = a + n*k
  with a= opening penalty, n= extension penalty, k= length of the gap

Compared to the linear gap penalty, this will give a small bonus to longer consecutive gaps compared to skipping one nucleotide, matching one, skipping one again etc (because in the case of longer gaps, you only pay the opening penalty once, so the global cost per base decreases).

Cells on your diagram:
From what I understand, you're working with an affine penalty of g(k) = 1 + 1*k. For the [3,4] you're "coming" from the score of 2 on the left (cell 1,2). You either:

match the A (row 1) with G (column 3) and get a penalty of 1 for substitution -> total of previous 2 + new 1 = 3
open a gap to skip the G (column 3) and get a penalty of 2 (1 for opening + 1 because gap length of 1) -> total of previous 2 + new 2 = 4.

So you choose the lesser penalty and go with 3, which is the number selected in the bigger box.
The explanation is the exact same one for the other cell you mentioned.


Remapping genomic coordinates to account for indels

I think that you need a LiftOver Chain file to transform your coordinates. You can obtain such a file using bcftools consensus with the -c parameter:

-c, --chain <file>         write a chain file for liftover

Then you can use it to transform coordinates in various genomic formats using CrossMap.


Visualizing protein-ligand interactions with LIGPLOT or any other suitable software?

You can visualize with Discovery Studio Visualizer and Schrodinger.
Arpeggio and PLIP can generate output files of the different interactions.


How to introduce point mutations in a protein structure?

PDBFixer and pyMutate are simple, lightweight, free and open-source. Both don't do any optimizing, though, so you may have to run some molecular dynamics afterward to get what would be the right "orientation".
VMD would work for sure, I believe PyMOL too. They are both graphical interfaces but you can start them in text only mode (or just import pymol in a python script and use the Python API). They are free for particulars and academics, but given their renown you may already be using them.
If you're looking for something a tad more complete and that may adapt the side-chain to the surrounding environment I'd go with Modeller, which provides two example scripts to do point mutations (here and there).
If you don't want to just get the side-chain of this amino-acid right but actually strive to see some general backbone adaptation too then use PyRosetta. It can probably do what all the others do, too, but given it is 1 GB big it may be overkill. Any of the package coupled with molecular dynamics would also do the trick.


How do I pull singe cell RNA sequencing data from GEO database?

The simplest would be using a count matrix (at the end of the link you have shared, section "Supplementary file"). For example GSE138651_Vagal_WholeNodoseSeq_raw.mtx.gz corresponds to "raw" count matrices obtained after mapping and counting the reads. This file, along with the GSE138651_barcodes.tsv.gz and GSE138651_genes.tsv.gz can be used to "read" the data into R; to that end there are many options in the form of packages, Seurat is one of the best documented and hence is beginner-friendly.


Manually define clusters in Seurat and determine marker genes

I think you are looking to FindAllMarkers function from Seurat. As you said, you just have to define your ident, that have to have the structure of a table (cell names as names and cluster as value):

pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

And then run the FindAllMarkers function:

FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)

You can specify several parameters in this function (type of DE to perform, thresholds of expression, etc).


Chromatin accessibility level on promoters with different CpG ratio calculation

Below I import the peaks call in narrowbed format from macs2 into a GRanges object. If you have a bed file, you can just use rtracklayer for it:
library(rtracklayer)
library(BSgenome.Hsapiens.UCSC.hg38)

peaks_narrowbed = "https://www.encodeproject.org/files/ENCFF846JAO/@@download/ENCFF846JAO.bed.gz"

extraCols <- c(signalValue = "numeric", pValue = "numeric",
                          qValue = "numeric", peak = "integer")
gr <- import(peaks_narrowbed, format = "BED",extraCols = extraCols)

head(gr)

GRanges object with 6 ranges and 6 metadata columns:
      seqnames        ranges strand |
         <Rle>     <IRanges>  <Rle> |
  [1]     chr1   10049-10326      * |
  [2]     chr1 180686-181040      * |
  [3]     chr1 186648-186903      * |
  [4]     chr1 191293-191604      * |
  [5]     chr1 267872-268146      * |
  [6]     chr1 586110-586324      * |
                                                   name     score signalValue
                                            <character> <numeric>   <numeric>
  [1] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_1       135     7.10907
  [2] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_2       309    11.65888
  [3] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_3        81     5.40289
  [4] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_4        56     4.54981
  [5] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_5       145     7.39343
  [6] atacseq_170613.t8_rep1.dedup.masked.sorted_peak_6       116     6.54035

Then we extract the sequences that are within these regions, you can of course redefine the region by extending it from the peak:
library(Biostrings)
SEQ = getSeq(Hsapiens,gr)
head(SEQ)

A DNAStringSet instance of length 6
    width seq
[1]   278 TAACCCTAACCCTAACCCTAACCCTAACCCTAAC...CCCCAACCCCAACCCCAACCCCAACCCCAACCC
[2]   355 CTAACCCCTAACCCCTAACCCTAACCCTACCCTA...CTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGA
[3]   256 GCAACCACCTGAGCGCGGGCATCCTGTGTGCAGA...CCCCTTCTTTCCATTGGTTTAATTAGGAACGGG
[4]   312 CAATCAGCAGGGACCGTGCACTCTCTTGGAGCCA...CCCTGCCCCTGTCTCTTCCGTGCAGGAGGAGCA
[5]   275 TGGATGTGTGCATTTTCCTGAGAGGAAAGCTTTC...TTGTTTAGTTTGCGTTGTGTTTCTCCAACTTTG
[6]   215 GAGGCTGAAGGAGACTGATGTGGTTTCTCCTCAG...CCAACTTTGTGCCTCATCAGGAAAAGCTTTGGA

This is the CG, and you can decide how you want to bin your regions:
freq = dinucleotideFrequency(SEQ,as.prob=TRUE)
hist(freq[,"CG"],br=100)


Delly finds identical deletions with different consensus sequences

It turns out that this was just a bug which has been fixed by this commit. The current master version of delly no longer exhibits this behavior.

If you are facing this issue, either build the current master version or wait a bit for the next release which should also include this fix.


Introduce errors in reference transcripts according to external dataset error model

Do any of the answers for this question help? Karel Brinda has mentioned a few read simulators in the answer to that question, and has a thesis with more information.

Excluding INDEL errors doesn't sound like a good idea; length can still be preserved even if doing that, it just needs an adjustment at the end of the sequence. Note that if you're trying to model nanopore reads, what you're really modelling is the base-caller, rather than the sequencer. I mention this in more detail in my answer.

In most cases where errors are modelled, I find it better to use publicly-available data instead. Especially for nanopore data, there are unmodelled systematic errors in the base-callers and sequencer that can't be simulated using any programs (because they are unmodelled). The following paper would be a good place to start for cDNA sequences, which looks at single-cell data from mouse (C57Bl/6) B1a cells:

http://www.biorxiv.org/content/early/2017/04/13/126847

Illumina and ONT reads for that study can be found in SRA under accession number SRP082530.

I don't know of any recent D. melanogaster studies that have been done using nanopore. There's always the option of spending $1000 on a purchase of a MinION with an RNA starter kit to do the study yourself. Here's an older targeted gene study, but bear in mind that it was using an R7.3 flow cell, so errors rates will be much higher than what is currently available:

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0777-z


One possible workaround would be to blast your sequences against the genome. Then, when you have found where they align, you can extract those regions only into a new fasta file (make sure to take a few kb around each target region) and then use that fasta file as the "genome" you pass to bowtie2.

That would give bowtie a far smaller file to work with and should reduce the memory needed significantly.


How exactly is "effective length" used in FPKM calculated?

The effective length is $\tilde{l}_i = l_i - \mu + 1$ (note the R code at the bottom of Harold's blog post), which in the case of $\mu < l_i$ should be 1. Ideally, you'd use the mean fragment length mapped to the particular feature, rather than a global $\mu$, but that's a lot more work for likely 0 benefit.

Regarding choosing a particular transcript, ideally one would use a method like salmon or kallisto (or RSEM if you have time to kill). Otherwise, your options are (A) choose the major isoform (if it's known in your tissue and condition) or (B) use a "union gene model" (sum the non-redundant exon lengths) or (C) take the median transcript length. None of those three options make much of a difference if you're comparing between samples, though they're all inferior to a salmon/kallisto/etc. metric.

Why are salmon et al. better methods? They don't use arbitrary metrics that will be the same across samples to determine the feature length. Instead, they use expectation maximization (or similarish, since at least salmon doesn't actually use EM) to quantify individual isoform usage. The effective gene length in a sample is then the average of the transcript lengths after weighting for their relative expression (yes, one should remove $\mu$ in there). This can then vary between samples, which is quite useful if you have isoform switching between samples/groups in such a way that methods A-C above would miss (think of cases where the switch is to a smaller transcript with higher coverage over it...resulting in the coverage/length in methods A-C to be tamped down).


Error Code in WGCNA: Dealing with large data sets: block-wise network construction andmodule detection

The problem is that blockwiseModules split your data into (probably 3) blocks, because the default maxBlockSize is 5000 and smaller than the number of genes in your data. If you have enough RAM (16GB is recommended for 12k genes, but 8GB just might be enough), rerun blockswiseModules again with maxBlockSize raised to above the number of your genes. Since you were able to run the function with maxBlockSize=8000, hoepfully you'll be able to do the single block calculation.


What is a signal 'chunk' in the context of Nanopore sequencing?

If you have access to the nanopore community website you can read more about this here: https://community.nanoporetech.com/posts/not-chunky-enough (I'm not sure to which extent I can copy paste from there so I won't).

A TLDR version: basecalling doesn't happen simultaneously for an entire read, but it is broken in segments which are processed separately and combined afterward. Those segments are chunks. Basecalling results are affected mostly by low-complexity regions, in which a heavily biased region (in terms of nucleotide content) leads to artifacts in the per-chunk basecalling. Those issues are solved by taking a chunking-size larger than most low-complexity regions.


Not exactly what you were looking for but may be you can use something like this?
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" \
    | xtract -pattern INSDSeq_feature-table \
        -group INSDFeature \
        -if INSDFeature_key -equals 'mat_peptide' \
        -tab '\n' \
        -element INSDFeature_location,INSDQualifier_value \
    | grep 'Segment 2'

For explorations of the INSDSet XMLs, xtract has -insd option that you can use. In this particular case, you can do something like this:
efetch -db nuccore -id "$ncbi2" -format gpc \
    | xtract -insd mat_peptide product peptide

For more information, checkout EDirect documentation and some  example queries


umap and Louvain clustering on normalized data

If the data are single-cell RNA seq, after normalization can we do
  UMAP and Louvain clustering ?

Typically people run PCA, UMAP and Louvain clustering on the normalised and log-transformed expression counts (but do marker gene and differential expression analysis on the non-normalised values). However, this remains controversial. I recommend reading "Current best practices in single‐cell RNA‐seq analysis: a tutorial"- it's a bit outdated, but still a good reference.

  By adequate I mean the clusters are the same but some are split into
  two, which makes sens looking at other results (like transcription
  factor analysis).

The number of (or size of) clusters should not be used to interpret how good the clustering methodology is, since it is generally dependant on the input parameters, e.g. a higher resolution parameter leads to more clusters. Instead, display the clusters on the UMAP and see if they are continuous 'visually' and check what genes are differentially expressed between clusters and see if they make sense 'biologically'.


bed file merge book-end features only if same name in column 4

awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed

That produces:

chr1    0   10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   31400   MIR1302-2HG::chr1:29553-31109.E1

This explicitly checks for bookends and nothing else and assumes that everything is sorted.


What is the best way to map/align my reads on a given genome?

This is a pretty straightforward analysis, assuming that you want to measure something like "how many reads align to genome 1, and how many reads align to genome 2?".

The bowtie2 manual has a section on this.
There are other tutorials for bowtie2.
bwa is also well-suited to this task.

You should be able to then use samtools to compute statistics on your alignments.


If you print G_list, you'll see this half way down:

1 ENSMUSG00000060550      15018      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]
2 ENSMUSG00000060550     110558      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]

So there are two entrez ids that map to the same Ensembl id. One of them maps to an alt locus.

If you exclude "entrezgene" from the query, you'll get 30 results.


Can pair end reads with high heterozygosity be used to polish PacBio assembly?

Resequencing your sample is not going to make it any less heterozygous. I am assuming that you are working with a diploid organism? Please update your question and specify the genome size and ploidy of your organism (the type of organism will help us update our answers, too).

The answer, yes, you can use Pilon + short reads to correct the PacBio assembly even with a highly heterozygous sample. To do so you should use the --diploid option in Pilon to let it know that there will be heterozygous sites. The way that it corrects heterozygous SNPs and heterozygous indels is different, and it seems to perform best with SNPs and sometimes does not correctly identify heterozygous indels (source).

Pilon itself doesn't handle phasing, but just corrects with whatever reads are in the bam file that it is provided with. The end result of providing a bam with reads from both haplotypes is that the corrected sequence will be a mosaic of the dominant allele at the heterozygous site as was observed in the short read data. This may or may not select for some alleles based on sequencing coverage bias of Illumina sequencers.


how do i change a variable with iterator under a for loop in python

cluster_0 to cluster_12 are just strings. they cannot be the names of arrays or matrices. Strings can only be used in print statements and string processing statements.

The correct way to define an array or matrix is as follows: call the arrays as: cluster[0] to cluster[12]

then, assign:

  for i in range(12):
        cluster[i] = gene_list.gene[gene_list.cluster==i]

then whenever you need to print the names of these arrays, you use:  "cluster_"+str(i)


Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV.

Update,
I can provide an answer. It's HDF5 with pandas dataframes (Python).

Simply

df.to_hdf('path/file.h5', 'key', format='table')

This isn't SQL but it works:

What it gives,

it's fast, that's it's purpose (not my main reason though)
it perfectly preserves all pandas dataframe information, this is  important
The 'key' ( object) allows you to store multiple dataframes into the same file, you just provide different key names

The last point is very useful, because you can pull out a dataframe of choice and the SQL commands are mimicked within pandas

I ran into trouble using CSV to store/recover complex pandas dataframe formats.

So it's good enough and much faster than SQL. I could have problems with format changes in time (unlike SQL), but for the here and now it works well.

Just to add a caveat, for dictionaries the json dump works well, but it is just shunting dictionaries around, so isn't relevant.


Targeted NGS, up to 99% of reads have been marked as duplicates

If you're doing target capture, especially of a small region, it is entirely possible that there are only five unique reads. Target capture is molecularly finicky and sometimes the entire region is not captured.


You would need to do this when creating the bigWig file. Note that I'm not aware of anything that allows manually setting the various zoom levels (libBigWig/pyBigWig will allow you to set a max number, but that's it).

This then raises the question of whether IGV is using an appropriate zoom level to begin with. I've never looked at how it chooses which of the zoom levels to use, but I know it uses them since if you don't have any it will fail to load bigWig files.

Honestly, though, if it's taking a minute to jump to a new region then I suspect that the network is just slow or the remote server is overtaxed. Changing the internal indexing isn't going to fix that.


How to extract metadata from NCBI's short read archive (SRA) for a few runs?

Entrez Direct is perfect for this sort of thing. You can do something like this:

$ esearch -db sra -q 'SRR390728' \
  | efetch -format runinfo \
  | cut -f1,12 -d,
Run,LibraryName
SRR292241,HS0798
SRR390728,HS0798

The runinfo table format has 47 fields including experiment accession, project accessions, bioproject and biosample identifiers and such. You may want to explore that a little bit more if you are interested in downloading metadata in addition to the LibraryName.

Also, if you prefer to parse XML output, you can append -mode xml to the efetch command to retrieve the runinfo table in XML format. xtract tool that comes with the Entrez Direct package can then be used to convert that into a table.


Biological replicates on Chip-seq Transcription factor data

You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only?  This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method.

When it comes to interpreting the biology behind the binding patterns, you are probably going to want to only interpret peaks which you find in both replicates.


WGCNA: Problem with selecting soft threshold

Don't worry about it too much. I would go with power 8 based on my general experience (also reflected in WGCNA FAQ) and on the mean connectivity around 50 and median around 20, which seems reasonable to me.


How does picard's MarkDuplicate handle unmapped reads?

From the Picard documentation:

DUPLICATION METRICS:
Metrics that are calculated during the process of marking duplicates within a stream of SAMRecords.

UNMAPPED_READS    The total number of unmapped reads examined. (Primary, non-supplemental)

It won't alter the flags on these reads, but it will count them in the summary report it generates.
You should be able to test this yourself with a small set of mapped + unmapped reads


Obtain CD (cluster of differentiation) marker expression information

I can offer this as a way to scrape the data from the Abcam resource:

library(rvest)
library(dplyr)

antigens <- read_html("http://www.abcam.com/primary-antibodies/human-cd-antigen-guide") %>%
  html_node("table") %>%
  html_table(header = TRUE)

# what targets expressed in "Langerhans" ?
antigens %>%
  filter(grepl("Langerhans", `Cellular Expression`)) %>%
  select(`Target Name`)

   Target Name
1         CD1a
2         CD1b
3         CD1c
4         CD1e
5         CD14
6         CD39
7         CD50
8         CD69
9         CD74
10        CD83
11       CD101
12       CD206
13       CD207

The data could be further processed by, for example, splitting the Cellular Expression column on either and or , (comma) to try and make one row per cell type.


Understanding DESeq2 design, contrast and results

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:

dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)

The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").

As an aside, this is probably a case where the edgeR-preferred way of creating groups of genotype-treatment combinations and then using a model of ~0+group might make your life a bit easier. You'll get the same results (more or less) regardless, but it'll probably be easier for you to think in those terms rather than remembering that the base level will be treatment HS30 and geno prg1.


biopython-motifs , How can I create motif from sequences of different length?

Motifs are not designed to be variable length in general since PSSMs/PVMs don't support insertions/deletions (indels). Instead, I suggest profile HMMs, such as those offered by HMMER3, which do support indels. You can parse HMMER3 output with Biopython via Bio.SearchIO.HmmerIO


Timeout when downloading the ncbi nr blast database

Same here. I used axel to download it a week ago, but fail to get finished. Also, the downloaded taxa zipped file (60Mb) doesn't match md5sum. Apparently, they update the FTP files every 30 minutes or 1 hour.


How to obtain NCBI Nucleotide (Entrez) IDs for all E coli sequences?

Do you mean the NCBI taxonomy IDs? If so, I've personally always preferred the way UniProt displays NCBI taxonomy, and it's easier to download IDs in mass.
For example, if you find E.coli's entry: https://www.uniprot.org/taxonomy/562
You can then click on "all lower taxonomy nodes", which will show all E.coli strains. If you click 'Preview first 10' before downloading, you'll see that the first column is the NCBI taxonomy ID of each strain of E.coli.
For example, E.coli strain 0.1197 is reported to have the ID 1005443 from this preview page. Indeed, if we search this ID on NCBI's taxonomy page, we do get the correct entry: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1005443
If this is what you want, click 'download all' to get a nice tab-separated file of all the strains of E.coli and their NCBI taxonomy IDs. This file could then be read and processed before you send API requests to Entrez.
Hope this helps!


How to detect mismatch before mapping in RNA-Seq data

At each step in the SLiTing and pooling the barcodes added are not random, but are taken  from a predefined list of possible barcodes.

When the reads arrive, first you filter by reads that have >1 low quality bases in the UMI region (which is truely random).

Then you compare the sequences in the 3 cell-barcode regions against the list of allowed barcodes for each of those positions. If a read is not within an edit distance of one from one of those barcodes, it is discarded.

Finally the reads that pass these filters can be aligned.


I want to write python program to cut a DNA sequence at an EcoRI restriction site and print the two fragments after cutting

import Bio
from Bio.Seq import Seq
my_seq = Seq("aaaaaaaaaaatttttttttttgaattccccccccccggggggggggg")
from Bio.Restriction import *
print(EcoRI.search(my_seq))
print(EcoRI.catalyse(my_seq))

Actully the EcoRI cuts the DNA after G in restriction site GAATTC.


Dealing with indels when converting a position to chromStart and chromEnd for a bed file

You can convert to VCF and use BEDOPS vcf2bed --insertions, vcf2bed --deletions and vcf2bed --snvs to segregate the VCF into three subsets.

For convenience, BEDOPS conversion tools deal with index shifts, including vcf2bed. The resulting BED file is 0-indexed, half-open.


How to identify gene expression signatures from gene expression data?

There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a gene signature:

Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down regulated genes and defines them as the signature of the group of cells. That's how the gene sets of the MSigDB which are not from (or derived from) pathway databases or gene ontologies are defined. There is no gold threshold of a fold change or a specific number of genes, it (highly) depends on the experiment you run and how good are the isolation of the cells.
Tools: DESeq2, limma, edgeR,...

As suggested by gringer, another way is to find which genes are expressed together (independently of how do they relate to other groups). This method relays on having a good enough number of samples from the same condition to be able to identify which genes are expressed the same way in those cells. If these genes are not expressed that way in other cell lines, then you have your signature for this group. Again, there isn't a golden method, different types of co-expression networks can be build and depends on how specific do you want to be or which is the data provided.
Tools: WGCNA, diffcoexp, TReNA,...

This is different from a gene set enrichment analysis1, which consists on evaluating if a group of genes (usually a gene set, a pathway or genes from a specific gene ontology) are not randomly sorted in a sorted list of genes (usually this sorted list is the fold-change of a differential analysis).

There are several methods developed depending on the null hypothesis, the type of sub-method and what do they take into account...
Tools: GSEA, PLAGUE, GSVA, limma, fgsea,...

TLDR: To perform a gene set enrichment you don't need a gene signature but a gene set (that can be a gene signature or a pathway or a random group).

1: Note that gene set enrichment analyisis is different from gene enrichment analysis. For a good review on the methods I recommend this article


How do I differentiate outliers from in-group variations in a DESeq2 PCA plot of 9 samples distributed into 3 conditions?

Since you have so few points here, it is hard to tell which are outliers by just looking at the PCA plot.
I suggest you use correlation heatmap for your gene expression and check whether there are obvious outliers.

As for using more than one groups when plotting PCA, I think you can just try passing a vector with many groups to the argument intgroup.


How to transform and sort the matrix to make a heatmap showing signatures?

You are almost there except for a few tricks:

If you want columns to be ordered using the values in Cell 4, you will need to provide a custom column order either with the column_order argument of the Heatmap() function or by providing ordered data. Below is the code to perform the latter. Most of the code is yours, I just tweaked it a little.

library(tibble)
library(dplyr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)

# Since you want to use cell4 as reference
cell4 <- as.data.frame(t(data["Cell4",]))

df <- cbind(cell4,
            Group = as.character(ann$Group))

# ordering by Group and by expression value, minus denotes decreasing order
df <- df[with(df, order(Group,-Cell4)),]

# to make it compatible to your code in your post
df <- rownames_to_column(df, "Samples")

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group" = c("A" = "orange", "B" = "darkgreen"))
colAnn <-
  HeatmapAnnotation(
    df = ann,
    which = "col",
    col = colours,
    annotation_width = unit(c(1, 4), "cm"),
    gap = unit(1, "mm")
  )

# rownames of the annotation df and the colnames of the data should be in
# the same order
data <- data[, df$Samples]

myCol <-
  colorRampPalette(c("#ffffff", "#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0, 1, length = 100)
hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = TRUE, # you asked your cells to be ordered
    cluster_columns = FALSE, # column order is the same as the provided data
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right"

)

The resulting heatmap, still not quite what you want:

You will realize that the values for Cell 4 are ordered from high to low for each of the two groups. To make Cell 4 be drawn at the top, you will need to reorder the rows. To do so, you will need to extract the row orders with the row_order() function:

> row_order(hmap)
 [1]  3  2 14  4 12 11  9 10  7  1  5 13  6  8

Now you can re-draw the heatmap after altering the row order somehow. For this sample case, doing it manually.

# promoting Cell 4 to the first position
modified_row_order <- c(4, 3, 2, 14, 12, 11,  9, 10,  7,  1,  5, 13,  6,  8)

hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = FALSE, # set to FALSE as a custom row order is supplied
    cluster_columns = FALSE,
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL,
    row_order = modified_row_order # argument to supply custom row order
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right")

If I have misunderstood your question and you do not want columns to be ordered according to expression values in Cell4, forget about the first part with the ordered annotation data frame.


How to provide a GFF/GTF file to Chromonomer

I think the part of the code that should allow you to do this is in chromomer.cc.

Lines 2857-2983 form the function parse_command_line. Nowhere is there any reference to anything resembling parsing a GFF file, so I think the documentation is not in line with the code at this point.

For reference, I looked at version 1.07. It would be helpful if the code were on github.


The Persistent uniform resource locator or PURL is one such solution, these are designed to be a bit more robust than permalinks in so much as they are supposed to survive the change of domain name.  The bio ontology community already use them http://purl.bioontology.org/docs/index.html


Differentially expressed genes analysis in Seurat

I'm not quite sure what you are asking, but you can use anything in the metadata as a grouping with FindMarkers.
DefaultAssay(combined.all) <- "RNA"
Idents(combined.all) <- "orig.ident"

for (sample in unique(Idents(combined.all))) {
  clustermarkers <- FindMarkers(combined.all, ident.1 = sample)
  print(paste("Cluster", sample))
  print(head(clustermarkers, n = 10))
}


Extract residue sequence from pdb file (in biopython, but open to recommendations)

There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

Try help(PDBParser) for more, or see:
http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#__init__


How to get transcriptome FASTA file for viruses for Kallisto pseudo-alignment?

You might be interested in the viral sequence data from CATCH. There are various different viral subsets available.
You could also try the NCBI virus website, filtering on human hosts. There's a "Download" link on this page that allows you to download the full sequence set:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?VirusLineage_ss=Viruses,%20taxid:10239&HostLineage_ss=humans,%20taxid:9605&SeqType_s=Nucleotide
You've asked a very specific question, so I'm giving you a very specific answer. It seems that what you actually want to do is different from the self-explanatory question.
For the purpose of creating an index for Kallisto or Salmon, it should be okay to use the full dataset. Those programs are designed to handle eukaryotic transcriptomes as input, which typically include many gene copies and isoforms, so transcript / genome duplication shouldn't be much of an issue.
Without knowing more about your particular application, it's difficult to provide additional advice. From what you've written, i.e. "what proportion of my raw data is originated from virus RNA", a hit to any viral sequence would be sufficient. For example, if you want any more specific information about viral proportion / species, then you should be using a metagenomic tool like Kraken2, rather than Kallisto / Salmon.


How to demultiplex pair-end fastq reads with barcode 2 in the identifier line?

You could try something like that:

#!/bin/sh

fq_in="${1}"

bioawk -c fastx '{print substr($comment, 7, 15)}' ${fq_in} | uniq \
    > "barcode1.txt"
bioawk -c fastx '{print substr($seq, 1, 5)}' ${fq_in} | uniq \
    > "barcode2.txt"

for bc1 in $(cat "barcode1.txt")
do
    for bc2 in $(cat "barcode2.txt")
    do
        fq_out="${bc1}_${bc2}.fq"
        bioawk -c fastx -v bc1=${bc1} -v bc2=${bc2} \
            'substr($comment, 7, 15) == bc1 && substr($seq, 1, 5) == bc2 {print "@"$name" "$comment"\n"$seq"\n+\n"$qual}' \
            ${fq_in} \
            > ${fq_out}
    done
done

This might be slow. You can a little bit of time if you directly provide lists of barcodes for the for loops.

It could be more efficient to write a custom program that can write to multiple files. It happens that I wrote something that does this in python just yesterday, but only reading the inline barcode. Modifying it to also use the in-header barcode is likely simple.


It is not like your conventional frequency which adds up to one. The density is low because the width of your bins is huge, and the number of observations you have is low. From this by Wickham, the basic kernel is:

where K is the kernel and h is the bandwidth (Scott, 1992b)

If you don't specify it, by default h will revert to nrd0:

library(rtracklayer)
gr_narrowPeak <- import("https://www.encodeproject.org/files/ENCFF072FBH/@@download/ENCFF072FBH.bed.gz",format = "BED",extraCols = extraCols_narrowPeak)
gr_chr21 = gr_narrowPeak[seqnames(gr_narrowPeak)=="chr21"]

da = as.data.frame(gr_chr21[rep(1:length(gr_chr21),gr_chr21$signalValue)])
dens = density(da$start)
dens$bw
[1] 1133889

plot(dens)

ggplot(da,aes(x=start))+geom_density()

You can decrease h such that you get something smaller (but at the cost of looking coarse). Important thing is it doesn't sum up to one.

If you would like to explain it easier like a histogram, consider converting it into counts, like explained in this post:

binW = 50000

ggplot(da,aes(x=start))+
geom_density(aes(y = ..density..*(nrow(da)*binW)),bw=50000) +
geom_histogram(binwidth=binW,fill="steelblue",alpha=0.5) +
theme_bw()

ggplot(da,aes(x=start))+
geom_histogram(aes(y=..count../nrow(da)),binwidth=50000) +
geom_density(aes(y = ..density..*binW),bw=binW,fill="steelblue",alpha=0.5) +
theme_bw()


Decontaminating RNA seqs for de novo transcriptome assembly and annotation of novel eukaryotes

according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them.
Try to run fastp with the following options:
fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

and run FastQc again to check to efficacity of the trimming.


Does the ".full.aln" file produced by snippy-core contain all bases of my input sequences aligned to the reference genome?

snippy-core is a tool that processes the snippy_outdir_{1..n} folders and produces the following files:

core.aln – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added
core.full.aln – a fasta format file that contains all the core sites and base calls at any site where the reference has a site. The reference is not padded but the sequences mapped to it may be padded with gaps and 'N's to match the length of the reference.  The alignment length does not change as more isolates are added.

snippy{-core} was written by Torsten Seemann et al.  snp-sites was written by Andrew Page et al.   snp-sites will reduce any alignment to the core polymorphic sites (columns) and will still output a fasta.  Therefore the output of a core.aln processed through snp-sites should be the identical the input.  All fasta alignments can then be processed downstream using your tree inference package of choice (e.g., FastTree, IQ-Tree, MrBayes, BEAST ...).  If you want a distance matrix output from a fasta alignment, use something like snp-dists, also written by Torsten Seemann, or the APE package in R.  Be sure to consider how gaps and Ns are treated during this step.


How to calculate average minor allele frequency (MAF) difference between two populations (i.e. European, Hispanic)

You need to report allele frequencies using plink for some European samples
Do the same for Hispanic samples
Use a software to subtract the two allele frequencies obtained above and obtain the result.


How to get the intersection of peaks after peak calling using MACS2?

bedtools intersect with multiple files after -b performs pairwise intersections between -a and each file in -b.
It does not intersect all files simultaneously.
See here, from the Bedtools documentation. The green section only contains intervals from -a that also intersect at least one interval in any -b file.

What I think you're looking for is bedtools multiinter (documentation here is unfortunately blank).
bedtools multiinter -i input [input...] will produce a BED file that covers every position listed in your input files, followed by how many and which of your input samples overlap that given position.

You can use this BED file and create something like an UpSet plot to visualize how many bases/peaks are shared between any combination of your samples.


How can I easily get the read size distribution of reads mapping on a certain set of regions?

This can be done rather simply by combining this samtools based answer and bioawk:

samtools view -L regions.bed aligned_reads.bam \
    | bioawk -c sam '{hist[length($seq)]++} END {for (l in hist) print l, hist[l]}' \
    | sort -n -k1

The samtools command extracts the alignements in the region given with the -L option, in sam format. The -c sam option of bioawk makes the $seq variable (corresponding to the read sequence field) available in the awk command.

The final sort ensures the size classes are ordered numerically.

We obtain a two-column size    count table.

Examples with timings:

$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | bioawk -c sam '{histo[length($seq)]++} END {for (l in histo) print l, histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.592s
user    0m0.816s
sys     0m0.036s

Since sam format is actually simple to parse, one can try to use mawk instead for slightly higher speed (the sequence being in the 10th field):

samtools view -L regions.bed aligned_reads.bam \
    | mawk '{hist[length($10)]++} END {for (l in hist) print l"\t"hist[l]}' \
    | sort -n -k1

(Contrary to bioawk, the default output field separator is not a tabulation.)

Example:

$ time samtools view -L /Genomes/C_elegans/Caenorhabditis_elegans/Ensembl/WBcel235/Annotation/Genes/protein_coding.bed results_Ribo-seq_test_1/bowtie2/mapped_C_elegans/WT_1/18-30_on_C_elegans_sorted.bam | mawk '{histo[length($10)]++} END {for (l in histo) print l"\t"histo[l]}' | sort -n -k1
18  201
19  560
20  1017
21  2568
22  5554
23  11221
24  27637
25  38797
26  17042
27  9853
28  15181
29  5888
30  528

real    0m0.362s
user    0m0.500s
sys     0m0.004s


Understanding the significance of BLAT score calculations for read fragment alignment

It would appear that BLAT and exonerate (as suggested by @terdon) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to BioPython's documentation, anyway:

  Most search programs only have HSPs with one HSPFragment in them,
  making these two objects inseparable. However, there are programs
  (e.g. BLAT and Exonerate) which may have more than one HSPFragment
  objects in any given HSP. If you are not using these programs, you can
  safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with BLAT and generating scores afterwards with the script described in their FAQ, or using exonerate which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!


Bioinformatics approach to dentifing potential PCR primer sequences for transcribed gene

It depends if you want to amplify your target from genomic DNA or cDNA.

Depending on intron size, primers that anneal to different exons may not be able to amplify across the intron(s), and if you're unlucky enough to design a primer that spans a splice junction, you'll definitely be stuffed if you're working with genomic DNA. I would try blastn or translate to an amino acid sequence and use tblastn to see if you can identify the intron/exon structure in a related species. If you don't have a genome sequence available for anything close, I think your only option is to target as small an amplicon as possible and hope for the best.

If you are working with cDNA, any standard primer design tool should do


How to efficiently compute the exact percentage of non-unique k-mers in a genome for given k?

You can do all of that with khmer. For example, abundance-dist-single.py produces a file with columns: k-mer abundance, k-mer count, cumulative count, and fraction of total distinct k-mers. So for question 1 you would sum column 2. For question 2 you would just get thek-mer countassociated with ak-mer abundance` of 1.

That package also provides a python API (see the read the docs link above) if you need to customize things more.


How to query a network for a bacterium, specifically Streptomyces caatingaensis or Streptomyces thioluteus?

I realized that some tools work better with Gene Names instead of some type of ID.

Using IIS I was capable of creating a network using this list of proteins (Streptomyces caatingaensis) by querying proteins by Genes Names having as base the E. Coli.

The tool has 3 main tabs: proteins, annotation, and interactome.
Search the Gene Names in proteins section, go to Interactome, set the parameters and create the network.

P.S. If you are not used to the tool you should read the research article available on the home page before using the network for research purposes.


It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing.

How to call variants (a vcf file) from mapped reads (a bam file) is very broad question, it depends on the sequencing technology (Illumina, PacBio, ...), the type of sequencing (whole genome sequencing, RAD, RNA-seq, ...), the sampling (calling of individuals vs population variability), the quality of a reference (human vs denovo assembly made by a postdoc in your lab)... For a lot of cases the answer is simply GATK, partially because it has well written manual, but mainly just because many people use it.

Some more detailed cases are well answered in more specific questions :

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
Somatic tumor only variant calling?
Why are there missing calls in a VCF file from exome sequencing?


Label protein strucuture residues in Chimera

You cant do this in Chimera, but in pymol you can write a script to do this.

In terms of the view, you can use get_view on the command line. This returns a section that you can copy and paste in - which will give you the same view every time. I recommend labelling the residue you want, then finding the position you want to view it from. Then using get_view to get that view, which you can use every time (just copy it into a text file) for a reference

This does presume your different molecules are aligned, but you can do that in pymol as well to make sure that when you use the output of get_view it is consistent.


How to use pyopenms.FeatureFindingMetabo() for peak picking in pyOpenMS?

This is a draft answer, if the OP still need needs I will extend it later.
Me too I found the documentation and examples of pyopenms very scarce. At the same time  OpenMS is huge, sometimes you find key information not in various other sources instead of the official documentation. This is especially true if it comes to metabolomics.
We built a Python OpenMS lipidomics LC MS/MS preprocessing pipeline, contained in this file:
https://github.com/saezlab/lipyd/blob/master/src/lipyd/msproc.py
The top class is the one called MSPreprocess, it implements the workflow, the other classes implement its steps.
Briefly, you start from profile mode mzML files, then

Peak picking with PeakPickerHiRes results centroided data
FeatureFindingMetabo assembles the traces in a metabolomics
compatible way
A MapAlignment algorithm alignes the traces across
scans
A FeatureGrouping algorithm alignes the features across multiple experiments

Our software mentioned above (lipyd) outputs data frames from the m/z, RT and intensity.
About OpenMS I can recommend to read the C++ code and docs. The Python classes and methods are automatically built minimal wrappers, everything works the same way as in C++.


What is the difference between NGS, 2GS, SBS and HTS?

Here are my attempts at definitions:

Sanger: A method of sequencing that depends on chain-terminatiing dideoxynucleotides. This sequencing uses the differential flow of DNA sequences of different lengths through a gel to determine the original DNA sequence, producing a single sequence per reaction container.

NGS: Next-generation sequencing, also referred to as 2GS (second-generation sequencing). This term is used to describe the first wave of sequencing technologies that followed Sanger sequencing technology. The use of NGS has become more confusing with the advent of long-read sequencing, because it's a common assumption that "next-generation" refers to the most recent technology (which is incorrect in this case).

HTS: High-throughput sequencing. This term describes any type of sequencing technology that produces large amounts of data, usually in the form of millions of different sequences produced from the same sequencing run.

SBS: Sequencing by synthesis. This term describes a method of sequencing that depends on the synthesis of [DNA] bases in order for sequencing to be carried out. This definition can extend into long-read sequencing (e.g. PacBio sequencers depend on synthesis during sequencing), but is more typically associated with only the second-generation sequencing technology.


The simple probabilistic scoring scheme you describe to query a target sequence against a (protein) family represented as a multiple sequence alignment (MSA) is actually similar in essence to a Position Weight Matrix (PWM), although PWMs are more refined. PWMs are powerful rapid tools often used for motif analyses.

However, your approach (and PWMs) don't handle insertions or deletions (indels), and you can't easily add these into your model (for example, consider that 2 1bp indels should not score the same as 1 2bp indel, and indels could occur within your MSA as well as in your target sequence). You could try and add some arbitrary gap opening and extension penalties somehow (like BLAST) but then you'd abandon your probabilistic model! By contrast, pHMMs have a fully fledged probabilistic indel model with insertion and deletion states.

Finally, of practical note, there are some excellent implementations of pHMMs for protein sequence analysis, such as HMMER3 and HHSuite.


why in 3' quantseq should we use reference genome and not transcriptome?

It would be great if you can elaborate on where you "found" this information.
In the transcriptome annotation you have a defined number of transcript per genes, but there might be others that are present in your sample and not in the annotation. For example you have treated the cells with some extreme condition that changes the 3' processing machinery.
With 3' quantseq, you are left with the read fragment that is at the end of your transcript. If the ends of these transcripts are not part of the sequences in the annotated transcriptome, these will be unaligned.
These instances might be rare if your transcriptome annotation is good, but let's say if you find that a gene of interest which should be expressed is not, and you suspect it might be due to a different isoform used, you can always check it with the genome alignment (see where reads might be falling in the gene). With a transcriptome alignment, this will not be possible.


Annotating conserved protein domains in bacterial assemblies

The first option. Using a pre-existing tool, such as PROKKA etc.

Diversity. Bacterial diversity is staggering and species such as E. coli and Salmonella, which are closely related diverged when mammals and dinosaurs/birds diverged, so searching nucleotide sequences is a bad idea as you will get nothing. Especially as you are hoping to find distant sequences that may be missed by the commonly used annotation tools.

Operons. Bacteria, especially thermophiles, often group their protein functionally in operons. Therefore, you can use the principle of guilt by association to guess what a mysterious protein does. It is extremely powerful. Say you have Voet & Voet handy, a pathway with a hole and a mysterious protein you can guess easily what it does.

Note, if you have a protein that you are particularly interested in, I would suggest making your own PSI PSSM file. The automated system cannot read the literature or do some quick checks to see when to stop, so may be overly conservative with the clusters.


Comparing AUCs: Discrimination of same Control from different Test Group - paired or unpaired? miRNA

Generally its paired. The bootstrapping personally I find concerning. I presume you are using it to augment your data set.
Bootstrapping at best tends to be a controversial statistical operation and could result in the non-significance you are observing. Anyway if its not significant between high and low - thats your answer.
You could of course perform a single ML classifying low, and high cancer risk against control. That would provide a method to understand the associations of miRNA.

Bootstrapping issue, I had a think about this. It depends on how the bootstrapping is done. In a paired sample this (I assume) samples with replacement at each paired position for both ROCs. If this is not the case - the test is junk. However, if this is the case, its actually a good test, because the paired structure of the data is always preserved no matter how many bootstrap samples are performed. So like if a 100 were done and the pairwise statistic (which we don't know what it is - but something like KS test) and its >0.05 for all 100 samples - I don't think there's any doubt thats a robust result.
Again if the two ROCs are bootstrapped independently, in my opinion the test is junk.
Reason The concern of bootstrapping is it causes weird distributions which were not present in the original data. This is only okay if there's lots of experience about interpreting the result. For example, in phylogenetics there's decades of experience interpreting bootstrapping and a consensus about what it means. This often isn't the case for most other areas of statistics.
I note this particular issue is really ramping up the rep score ;0)


Authoritative source on human cytogenetic regions?

Do you mean cytoband coordinates? For the human genome you can find them for example at UCSC.

hgdownload.cse.ucsc.edu/goldenPath/hg38/database/cytoBand.txt.gz


How to interpret amino acid representation

This means that at position 383 the Cysteine(C) is mutated to a Serine(S).  Depending on the organism and gene usually the mutation notation is linked to a specific transcript since the position could vary depending on transcript splicing.


Answer:
There are no methods that exist for investigating transcript data inside a cell without destroying the cell. Separating RNA from cells requires essentially destroying the cells and using a chemical cocktail to split its constituent components.

Alternative 1: One option is to investigate extracellular transcripts of living cells is to centrifuge cells in culture and isolating RNA from the growth media supernatant [1]. Cell-free DNA and cell-free RNA are hot topics in contemporary biology and bioinformatics as a potential source for disease biomarkers. [2]

Alternative 2: Another option to monitor a transcript from a living cell is to use a transgenic construct with a reporter molecule like luciferase or a fluorescent protein to produce light or become fluorescent when a promotor is activated [3, 4]. This does not directly measure transcript quantity but answers a similar question: how does transcription change in a living cell?

Sources:

[1] Théry, C., Amigorena, S., Raposo, G., & Clayton, A. (2006). Isolation and characterization of exosomes from cell culture supernatants and biological fluids. Current protocols in cell biology, 30(1), 3-22.

[2] Yu, X. M., Wu, Y. C., Liu, X., Huang, X. C., Hou, X. X., Wang, J. L., ... & Ling, Z. Q. (2016). Cell-Free RNA content in peripheral blood as potential biomarkers for detecting circulating tumor cells in non-small cell lung carcinoma. International journal of molecular sciences, 17(11), 1845.

[3] Kaskova, Z. M., Tsarkova, A. S., & Yampolsky, I. V. (2016). 1001 lights: luciferins, luciferases, their mechanisms of action and applications in chemical analysis, biology and medicine. Chemical Society Reviews, 45(21), 6048-6077.

[4] Rao, J., Dragulescu-Andrasi, A., & Yao, H. (2007). Fluorescence imaging in vivo: recent advances. Current opinion in biotechnology, 18(1), 17-25.


Find pattern that is present twice and allow <=2 mismatches on each

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

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

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

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

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

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


Random addition method for phylogenetic tree reconstruction

MEGA manual may help to explain the criteria used in selecting the final maximum parsimony tree.

The random tree additions refers to the generation of random initial trees for a heuristic search to find MP trees. In this case, a tree is generated by randomly selecting a sequence and adding it to the growing tree on a randomly-selected branch.

Instead of using the random tree additions you can use the min-mini algorithm for the initial tree set.


How to calculate Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation from protein (amino acid) sequences?

Rajarshi Guha describes the formulas to calculate Moreau Broto autocorrelation, Moran coefficient and Geary coefficient and has a bibliography.

But why reinvent the wheel? It seems that the Autocorrelation module of the PyBioMed package does already what you want: it is a pure Python implementation without any external dependencies.


How do I find a list of PDB structures bound to a peptide inhibitor or peptide-like inhibitor?

The quickest way to screen the PDB database is Advanced Search from either RCSB or PDBe.

Since you need entries with 2+ proteins you can set the number of protein chains to at least 2 in the RCSB search.  Or you can select "protein/protein complex" in Assembly composition in the PDBe search.

It's less clear how to search for inhibitors. If you'd look for the "inhibitor" word in the molecule names you'd get (together with the previous criterium) 1200+ entries.


We can skip the first row:

read.delim("txt.txt", skip = 1, header = FALSE)

Or we can get the "uncorrupted" version of the table, it looks like it is from this paper, Table S4.


You probably have red blood cell contamination.  It's probably better to identify which cells they are.  They'll probably cluster away from the cells you are interested in.


How to learn/get started bioinformatics?

The formal answer is O'Reilly "Bioinformatics programming, using Python". Check out the classic O'Reilly texts on "Mastering Perl for bioinformatics", and possibly"Beginning Perl for bioinformatics" (which is likely far too basic), which establish some core principles and have a solid reputation in teaching. However, the Perl/Python thing might get confusing (I dunno).

If you are statistics orientated at a guess, O'Reilly's, Bioinformatics data skills... I don't know, but something that would focus on R.

The real goal is core-technical competance because the biological understanding will be taught by your supervisor. Thereafter you might turn your attention to teaming up with a lab for project training.


Discordance in gene signature behavior between bulk and single-cell RNASeq

One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).

Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).

The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.


I have never used this program, but reading the manual in the link it has to mean the following.

The manual says:

  Writing -filt_rules_l 'eur.maf<0.05' on the command line would tell the program to remove any variants with eur.maf values less than 0.05 from the reference panel.

So, -filt_rules_l 'filter==0' would mean that the program removes any variants with filter values of 0 from the reference panel. In the legend file there should be a column named filter, which are filtered out when the value is 0.

Here is the example legend file, in your case there must be a custom additional column named "filter":

$ head ALL_1000G_phase1integrated_v3_chr1_impute.legend

id position a0 a1 afr.aaf amr.aaf asn.aaf eur.aaf afr.maf amr.maf asn.maf eur.maf
rs58108140 10583 G A 0.0407 0.1713 0.1311 0.2071 0.0407 0.1713 0.1311 0.2071
rs189107123 10611 C G 0.0142 0.0276 0.0140 0.0211 0.0142 0.0276 0.0140 0.0211


How to merge transcript sequence with same name in a FASTA file.

Certainly this python script is not the most elegant solution but it will give you the desired result:

import sys

seqs = {}
with open(sys.argv[1], "r") as fh:
    curr = ""
    for line in fh:
        if line.startswith(">"):
            curr = line
            if curr not in seqs:
                seqs[curr] = ""
        else:
            seqs[curr] += line.strip()

for ident, seq in seqs.items():
    print("{}{}".format(ident, seq))

If you saved the script in seq_merge.py just run it with:

python seq_merge.py your_file.fasta > result.fasta


Assessing PyMol sequence alignment object

In the application, unfortunately the only option is doing

cealign ,, object =
Menu > Export Alignment... on your cealign  object
and converting the Clustal format using one of the millions of tools online

get_fastastr
As you know cealign adds gaps, but cmd.get_fastastr('') strips them. There is no option to keep the gaps —annoyingly especially with missing loop info!
The help has this to say:
 get_fastastr(selection='all', state=-1, quiet=1, key='', _self=<module 'pymol.cmd' from '/Users/matteo/miniconda3/lib/python3.7/site-packages/pymol/cmd.py'>)

Correction and dissection

So I had a quick check to see if there is an alternative route to get the sequence string with gaps. These are my observations using the Python3 module, sorry if they aren't polished up.

In the comments I said something incorrect. cmd.cealign returns a dictionary with the 4x4 rototranslational matrix —no need for an object. Using Python3 and two TIM barrels as an example:
import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1TIM', 'tim')
    pymol.cmd.remove('not chain A')
    pymol.cmd.fetch('5AHE', 'hisA')
    x = pymol.cmd.cealign('tim', 'hisA')
    print(x)
# {'alignment_length': 184, 'RMSD': 3.8800425577361755, 'rotation_matrix': [...]}

While the object does not contain any useful info from PyMol directly*:
    ...
    pymol.cmd.cealign('tim', 'hisA', object ='rototrans')
    pymol.cmd.get_session('rototrans')['names'][0][5]
    # useful info is generally in ['names'][0][5]

(*) The code is actually wrapped C++ code, so the object rototrans is in the C++ layer and can only be access by session or possibly by some functions in cmd.exporting. In fact, it can be exported to an aln string (ClustalW).
    aln = pymol.cmd.exporting.get_alnstr('rototrans')

In the PyMOL application, doing the cealign with an object (cmd.cealign('tim', 'hisA', object ='rototrans'), allows the  Menu>Export Alignment...>'rototrans' which is saved as a ClustalW file, as you mention. Namely:
CLUSTAL

tim          APRKFFVGGN----WK-------MNGKRKSLGELIHTLDGAKLSADTEVVCGAPSI------Y
hisA         MIIPALDLIDGTVVRLHQGDYARQRDYGN---DPLPRLQDYAAQGAGVLHLVDLTGAKDPAKR
                           ..       ......   .....*..................      .

tim          LDFARQK---LDAKIGVAAQNCYKVPKGAFTGEISPAMIKDI-GAAWVILGHSERRHVFGESD
hisA         QIPLIKTLVAGVNVPVQVGGGVRTEE--------DVAALLKAGVARVVIGS---------TAV
             .......   ................        ..*..... .*..**..         ...

tim          ELIGQKVAHALAEG---LGVIACIGEKLDEREAGIT-----------EKVVFQETKAIADNVK
hisA         KSPDV-VKGWFERFGAQALVLALDV---RIDEHGTKQVAVSGWQENSGVSLEQLVETYLPVGL
             ..... *.......   ..*.*...   ...*.*..           .....*..........

tim          DWSKVVLAYEPVWAIGTGKTATP------QQAQEVHEKLRGWLKTHVSDAVAVQSRIIYGGSV
hisA         ----KHVLCTDI----------SRDGTLAGSNVSLYEEVCARYP---------QIAFQSSGGI
                 ........          .      .......*.......         *......*..

tim          TGGNCKELASQHDVDGFLVGGASLKPEFVDIINAKH------------
hisA         GDIDDIAALRGTGVRGVIVGRALLEGKFTVKEAIQ?WQNVKGHHHHHH
             .............*.*..**.*.*...*..

However, there is no method in cmd.exporting that gives out just the sequence with gaps.
Whereas the cmd.exporting has the following... none applicable.
 print(dir(pymol.cmd.exporting))
['DEFAULT_ERROR', 'DEFAULT_SUCCESS', 'QuietException', 'SDF', 'SDFRec', 'Shortcut', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_cmd', '_feedback', '_get_dump_str', '_get_mtl_obj', '_raising', '_resn_to_aa', '_session_convert_legacy', '_unit2px', 'absolute_import', 'assign_atom_types', 'cPickle', 'cache', 'cache_action_dict', 'cache_action_sc', 'cmd', 'copy', 'copy_image', 'export_coords', 'fb_mask', 'fb_module', 'get_alnstr', 'get_bytes', 'get_ccp4str', 'get_cifstr', 'get_fastastr', 'get_maestr', 'get_mol2str', 'get_pdbstr', 'get_pqrstr', 'get_psestr', 'get_sdfstr', 'get_session', 'get_str', 'get_xyzstr', 'io', 'is_dict', 'is_error', 'is_list', 'is_ok', 'is_tuple', 'loadable', 'lock', 'multifilenamegen', 'multifilesave', 'multisave', 'os', 'png', 'print_function', 'pymol', 're', 'save', 'savefunctions', 'selector', 'sys', 'thread', 'traceback', 'unlock']


Entrez.efetch returns incomplete genbank records

This has been answered before. See: https://stackoverflow.com/a/55402322/6262370

In short, you need to either use rettype='gbwithparts' or rettype='gb', style='withparts' to download the entire genbank flat file.


GSEA enrichr with 10x genomics differential_expression ranks

You should be using the most positive fold changes, not the most negative. The positive ones are the ones that are most highly expressed in your cluster of interest. You probably want to use less than 100 genes. Depending on the experiment, you may not have that many good markers for a particular cluster.
The other question you mentioned refers to DESeq2 specifically, so it is not applicable here.


How to convert fasta file to tab delimited file

If you have multi-line fasta files, as is very common, you can use these scripts1 to convert between fasta and tbl (sequence_name <TAB> sequence) format:

FastaToTbl

#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

TblToFasta

#! /usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

Save those in your $PATH, make them executable, and you can then do:

$ cat file.fa
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC
$ FastaToTbl file.fa
sequence1   ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2   GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

And, to get the Fasta back:

$ FastaToTbl file.fa | TblToFasta
>sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
>sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

This can be a very useful trick when searching a fasta file for a string:

TblToFasta file.fa | grep 'foo' | FastaToTbl

If you really want to keep the leading > of the header (which doesn't seem very useful), you could do something like this:

$ perl -0pe 's/\n//g; s/.>/\n>/g; s/$/\n/;' file.fa
>sequence1 ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGCTAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATAC
>sequence2 GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGCATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC

But that will read the entire file into memory. If that's an issue, add an empty line between each fasta record, and then use perl's paragraph mode to process each "paragraph" (sequence) at a time:

perl -pe  's/>/\n>/' file.fa | perl -00pe 's/\n//g; s/.>/\n>/g; s/$/\n/;'

1Credit to Josep Abril who wrote these scripts more than a decade ago.


gffread: GFaSeqGet errors on coordinate overhang

Basically, the source of the problem is issues with hard/soft clipped bases in short-reads RNAseq reads alignments to reference genome in post-transcriptome assembly steps.

Truncating the coordinates in the gtf file to match the reference is a solution. For this, I will share a perl solution (non-official) that will do the job.

perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf > transcripts_fixed.gtf

Where:

genome.fa.fai is the fasta index of the genome file
transcripts_to_fix.gtf is the transcript set out of the assembler with issues
transcripts_fixed.gtf is the new transcript file with the repaired end coordinates

Content of cufftrim.pl:

#!/usr/bin/perl
use strict;
use Getopt::Std;

my $usage = q/Usage:
 cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
umask 0002;
getopts('o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
if ($outfile) {
  open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
  select(OUTF);
  }
# --
my $fai=shift(@ARGV) || die ("$usage\nSamtools fasta index required!\n");
my $gtf=shift(@ARGV) || die ("$usage\nCufflinks' transcripts file required!\n");
my %slen;

open(FAI, $fai) || die("Error opening $fai\n");
open(GTF, $gtf) || die("Error opening $gtf\n");
while (<FAI>) {
 chomp;
 next if m/^#/;
 my @t=split(/\t/);
 die("Parsing error: $t[1] length for sequence $t[0] ?\n")
    unless $t[1]>10;
 $slen{$t[0]}=int($t[1]);
}
close(FAI);

# now process the GTF file
my ($tcount, $ecount); #trim counts
while (<GTF>) {
 my @t=split(/\t/);
 my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
 if ($t[4]>$l) {
   $t[4]=$l;
   $t[2] eq 'exon' ? ++$ecount : ++$tcount;
 }
 print join("\t", @t);
}
close(GTF);
print STDERR "$tcount transcripts, $ecount exons adjusted.\n";
# --
if ($outfile) {
 select(STDOUT);
 close(OUTF);
 }

Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution.


Why are my genes filtered for Gene Ontology term enrichment?

What is "click here for details" saying? That link should contain an explanation of what happened.

GSEA usually requires entrez_id to run. If your list contains gene_symbols, then before enrichment, you need to convert the genes names. This can and is (always) tricky as genes are known with different names. For example, by googling I couldn't find any AASS3 gene, but I found AASS, which is also named LKRSDH,LORSDH, and LKR/SDH. Depending on the software you are using for conversion, and the genome version used some genes will be missed. I usually run (in R) a combination of geneSynonym and biomaRt to fill the gaps in the lists.

Even better, would be to use ensembl_id or even better entrez_id directly


Common gene naming conventions and how to convert between them

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.


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

will index your .bam file.


How do I plot clusters based on DNA sequence alignment?

I think one of the reasons you struggle is that clustering DNA sequences is not a clearly defined task. In general intention of clustering is to reconstruct, or approximate the relatidness of the DNA seqeunces.
If all these sequences are homologous you can do a multiple sequence alignment (using MAFFT or Clustal) and build a phylogenetic tree (using something like IQ tree). If your prof suggested ClustalW, this is probably what they meant. Evolutionary speaking, that's the optimal method for getting relations of sequences, but they must be "alignable". Also with this many sequences, it will be computationally quite heavy.
For within species relatedness, phylogenetic reconstruction is sometimes impossible, because recombination breaks down the linkage of individual alleles. However, if all of the sequences are long enough, you can do a principal component analysis on SNPs, for example using something like plink. The plot you showed actually looks a lot like PCA (it even has the explained variability on axes).
Finally, regardless of their relationship, you can always do a dendrogram based on a similarity matrix, but for that, you need to define similarity. That will be usually something kmer-based. I am not following this field very closely lately, but there are some recent papers about kmer-based clusterings, like this one.
When looking for the paper I also run into this one, which proposes genome clustering using "Inter-nucleotide Covariance". No idea how useful/sound it is. Just found the abstract quite interesting.
Long story short, it really depends on what you are after and what are your sequences! If you can build a reliable tree, it is the most powerful way to determine the relationships among sequences.