# Question
Title: Question about the practical relevance of chromosome position and p-value
Article of interest: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3023908/#SD1

Referring to Figure 2b of that article, if I understand correctly, the x-axis refers to the position of each variant in a chromosome while the y-axis is a log transform of the raw p-value of association of each variant with some variable of interest.

My questions is, is there practical importance in finding strong association (high log transform p-value) around one particular position versus finding it on some other position in the same chromosome?

My follow up question is, are these log transformed p-values publicly available?

Context: I am a statistics student studying multiple testing. My biology background is very limited.

# Answer
In the specific case that you are talking about, the answer is very much yes that the specific position is important. Genes occur at (fairly) fixed positions along chromosomes.

The point of that figure, as indicated in the legend, is that the very low p-values are around the gene *RPM1*. That gene is known to have important natural genetic variation associated with the phenotype in question (the one for whichg they are looking for genetic associations). That figure is therefore an important validation of the methods of the authors.

In terms of finding such data, a fun account on twitter is the GWASbot, that just posts random manhattan plots (which is the informal name for the plot you are talking about), along with links to source data. I believe that that source data includes the p-values for each ascertained position, along with other summary statistics. One example is here.

If you are specifically interested in the topic of multiple test correction, an important related concept in GWAS is that of genome-wide significance.

**Update:**

To unzip and look at the file, use bgzip:

```
$ bgzip -d ~/Downloads/continuous-20491-both_sexes.tsv.bgz
$ head ~/Downloads/continuous-20491-both_sexes.tsv
# prints the following:
chr pos ref alt af_meta beta_meta   se_meta pval_meta   pval_heterogeneity  af_AFR  af_CSA  af_EUR  beta_AFR    beta_CSA    beta_EUR    se_AFR  se_CSA  se_EUR  pval_AFR    pval_CSA    pval_EUR    low_confidence_AFR  low_confidence_CSA  low_confidence_EUR
1   11063   T   G   2.088e-02   2.247e-02   2.072e-01   9.137e-01   0.000e+00   2.088e-02   NA  5.212e-05   2.247e-02   NA  1.287e-01   2.072e-01   NA  5.080e-01   9.137e-01   NA  8.000e-01   false   NA  true
1   13259   G   A   2.805e-04   9.840e-02   8.909e-02   2.694e-01   1.000e+00   NA  NA  2.805e-04   NA  NA  9.840e-02   NA  NA  8.909e-02   NA  NA  2.694e-01   NA  NA  false
1   17641   G   A   8.908e-04   -1.244e-02  4.922e-02   8.005e-01   1.000e+00   NA  NA  8.908e-04   NA  NA  -1.244e-02  NA  NA  4.922e-02   NA  NA  8.005e-01   NA  NA  false
1   30741   C   A   NA  NA  NA  NA  NA  9.400e-04   NA  NA  3.935e-03   NA  NA  1.297e+00   NA  NA  9.976e-01   NA  NA  true    NA  NA
1   51427   T   G   NA  NA  NA  NA  NA  1.127e-03   NA  NA  8.780e-01   NA  NA  9.366e-01   NA  NA  3.486e-01   NA  NA  true    NA  NA
1   57222   T   C   7.078e-04   -4.153e-03  5.656e-02   9.415e-01   1.000e+00   NA  NA  7.078e-04   NA  NA  -4.153e-03  NA  NA  5.656e-02   NA  NA  9.415e-01   NA  NA  false
1   58396   T   C   2.735e-04   3.011e-02   8.221e-02   7.142e-01   1.000e+00   2.622e-03   NA  2.735e-04   -1.069e-01  NA  3.011e-02   7.079e-01   NA  8.221e-02   8.800e-01   NA  7.142e-01   true    NA  false
1   62745   C   G   NA  NA  NA  NA  NA  2.701e-03   NA  NA  5.524e-01   NA  NA  5.529e-01   NA  NA  3.178e-01   NA  NA  true    NA  NA
1   63668   G   A   NA  NA  NA  NA  NA  NA  1.129e-03   3.673e-05   NA  2.205e-01   -4.462e-02  NA  9.738e-01   2.544e-01   NA  8.209e-01   8.608e-01   NA  true    true

```
> 1 votes

# Answer
> Are these log transformed p-values publicly available?

For human population data, it's pretty common for the association statistics to be either included with the research paper, or uploaded to a public database (e.g. see curated compilations like GWAS Central or NHGRI GWAS Catalog).

For plant and other non-human data, there doesn't seem to have been as much drive to include GWAS data with the research paper, and it's often up to the authors to decide to do that. I find this a little odd, given that there should be fewer ethical barriers to releasing non-human public data.

If you're interested in a particular dataset, you can contact the authors, but I acknowledge that such an approach can be too high of a barrier for most interested researchers.

> Is there practical importance in finding strong association (high log transform p-value) around one particular position versus finding it on some other position in the same chromosome?

It's common for statistical tests for association to assume independence of variant locations. With this in mind, two associated regions that are one base apart should not be any more interesting than two associated regions that are 50 megabases apart.

The problem with this assumption is that DNA variation is *not* statistically independent. Recombination of DNA (and associated reassortment of genetic variation) between two locations is more common when the locations are close, and this *linkage* should really be accounted for when doing association tests. One way to do this is to combine variants that share the same ancestral origin into a haplotype, then do association tests on the *haplotype block* rather than the single-base *genotype*. Even if the particular observed variants are not directly causative, there may be another untyped causative variant (or combination of variants) on the same haplotype block, meaning that the block becomes a proxy for the causative variant.

But this haplotype-based approach becomes quickly computationally intractable when looking at larger regions of chromosomes, and is very difficult to use for comparing populations with substantially different ancestral backgrounds. The workaround (if it could be called that) is to ignore the linkage aspect of genetic variation, and return to those same statistical tests, but declare regions that have peaks of close association to be statistically interesting, and zoom into them to look in more detail at the region.

That zoom in would be a great opportunity to use more complex / accurate haplotype-based tests for confirming an associative signal around a particular region. But for the most part, as with the paper you've indicated, the tests are identical.

However, even ignoring the haplotype issue... *\<begin rant\>*

The most common manhattan plots represent *only* p-value on the Y-axis. I think this is a mistake, because the p-value only represents the confidence in a statistic, and the important thing for association is the statistic itself, rather than the confidence in that statistic.

I have a few specific concerns about this pure p-value approach:

* there may be values that have less confidence, but a much higher associative statistic, and contribute more to the trait of interest (e.g. see here)
* values with a low association statistic are more sensitive to genotyping error, leading to spurious associations.
* the p-value doesn't indicate the direction of association
* p-values should not be used for ranking, and this ranking is implicitly done when using the p-value only manhattan plot

In addition, the generated p-values often don't make sense, and I haven't seen any commentary on this. Anything smaller than $10^{-40}$ doesn't represent any reasonable natural random occurrence, even on galactic timescales, but often these astronomically-small p-values can become insignificant through a few random mutations in one of the populations.

As a quick-fix using the existing statistics that are typically presented in GWAS results, where single points are desired, I recommend that people instead use the standard deviation to determine the minimum expected beta score:

see here for more details.

More rigorously, I think that people should consider doing their own internal validation of associations via bootstrap sub-sampling, and report a non-parametric statistic: the best observed *rank* after multiple bootstrap sub-samples. This requires access to genotype-level data, so cannot be done on existing GWAS results. I explore this in more detail in this preprint.

Yes, I'm aware this goes against the grain of how GWAS is typically done. I have been emboldened in my position on this after reading the ASA principles about the p-value, and seeing that many of the principles are in direct opposition to the typical manhattan plot:

https://www.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#\_i31
> 1 votes

---
Tags: statistics, gwas
---


# Question
Title: loading in .txt files with different spacing in column names into R
Trying to read in this file but saved as a .txt file: http://www.informatics.jax.org/downloads/reports/HOM\_MouseHumanSequence.rpt

`gene_conversion = read.table("HOM_MouseHumanSequence.rpt.txt",header=T, fill=T)`

but the column names are splitting weirdly.

```
> head(gene_conversion)
        DB  Class        Key Common Organism       Name        NCBI Taxon    ID                      Symbol
1 38710983 mouse, laboratory  10090     Gdnf      14573  MGI:107430 Chr15   3.8                          cM
2 38710983  human       9606   GDNF     2668  HGNC:4232 OMIM:600837  Chr5 p13.2   Chr5:37812677-37840044(-)
3 38710984 mouse, laboratory  10090     Rrm1      20133   MGI:98180  Chr7 54.72                          cM
4 38710984  human       9606   RRM1     6240 HGNC:10451 OMIM:180410 Chr11 p15.4    Chr11:4094685-4138925(+)
5 38710985 mouse, laboratory  10090   Hoxd12      15432   MGI:96204  Chr2 44.13                          cM
6 38710985  human       9606 HOXD12     3238  HGNC:5135 OMIM:142988  Chr2 q31.1 Chr2:176099802-176100760(+)
                                                            EntrezGene
1                                             Chr15:7840327-7867056(+)
2 NM_199231,NM_001278098,NM_199234,NM_001190469,NM_001190468,NM_000514
3                                          Chr7:102090902-102118978(+)
4                     NM_001033,NM_001318065,NM_001330193,NM_001318064
5                                            Chr2:74505357-74508049(+)
6                                                            NM_021193
                                                                                ID.1
1                                     NM_001301357,NM_010275,NM_001301333,NM_001301332
2 NP_001265027,NP_000505,NP_001177397,NP_001177398,XP_016864826,XP_011512332,NP_954701
3                                                                            NM_009103
4                                     NP_001024,NP_001317122,NP_001304994,NP_001304993
5                                                                            NM_008274
6                                                                            NP_067016
                                             Mouse    MGI ID.2 HGNC ID.3 OMIM Gene ID.4 Genetic Location Genomic
1 NP_001288286,NP_001288262,NP_034405,NP_001288261 P48540   NA   NA   NA   NA   NA   NA      NA       NA      NA
2                                           P39905          NA   NA   NA   NA   NA   NA      NA       NA      NA
3                                        NP_033129 P07742   NA   NA   NA   NA   NA   NA      NA       NA      NA
4                                           P23921          NA   NA   NA   NA   NA   NA      NA       NA      NA
5                                        NP_032300 P23812   NA   NA   NA   NA   NA   NA      NA       NA      NA
6                                           P35452          NA   NA   NA   NA   NA   NA      NA       NA      NA
  Coordinates X.mouse. X. human. X..1 Nucleotide RefSeq IDs Protein RefSeq.1 IDs.1 SWISS_PROT IDs.2
1          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
2          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
3          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
4          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
5          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA
6          NA       NA NA     NA   NA         NA     NA  NA      NA       NA    NA         NA    NA

```

# Answer
> 0 votes
I see spaces in the headings, which suggests this might be a tab-separated file, rather than a space-separated file. If that's the case, you should be using `sep="\t"` in your `read.table` command, or alternatively using `read.delim`:

```
gene_conversion <-
  read.table("HOM_MouseHumanSequence.rpt.txt", sep="\t", header=T, fill=T)

```

or:

```
gene_conversion <- read.delim("HOM_MouseHumanSequence.rpt.txt")

```

---
Tags: r
---


# Question
Title: calculating kmer nucleotide frequency per column
I have a list of sequences:

```
CAGGTAGCC
CCGGTCAGA
AGGGTTTGA
TTGGTGAGG
CAAGTATGA
ACTGTATGC
CTGGTAACC

```

Each sequence is nine nucleotides long. I want to calculate the dinucleotide frequency of AA, CG ..etc for each of the dinucleotide positions. Similar problem shared with one nucleotide in this post.

For example, 1st two column will have these bases CA,CC,AG,TT etc, and second two should be AG,CG,GG,TG and so on. Output should be something like this

```
AA   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
CC   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
GC   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
TC   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36

```

AA, CC, GC and TC could be any dinucleotide combination.

**What I tried from the last post: mainly taken from @Devon**

```
import itertools
hl = []

#making all possible dinucleotide combinations
bases=['A','T','G','C']
k = 2

kmers = [''.join(p) for p in itertools.product(bases, repeat=k)]
kmers = dict.fromkeys(kmers, 0)
print(kmers)
for i in range(9):
    hl.append(kmers)

# my seq file
f = open("H3K27ac.seq")
nLines = 0

for line in f:
    id = 0
# trying to read dinucleotide at each loop
    for (fst, sec) in zip(line[0::2], line[1::2]):
        id = id + 1
        hl[id][fst+sec] += 1
    nLines += 1
f.close()

nLines = float(nLines)
for char in ['AA', 'TT', 'CC', 'GG', 'CG', 'GC']:
    print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))

```

By somehow it's estimation the total count of all dinucleotide in whole file. I am trying to make it more general so that it can be used for tri- and any further number of nucleotides.

# Answer
> 2 votes
This should do the trick. The following code works for any aligned character set, but assumes that all sequences should be the same length.

```
import itertools

with open("f.txt") as seqs:
    sequences = list(line.strip() for line in seqs.readlines())

# Just in case sequence lengths are different
n = min(len(s) for s in sequences)
if any(n != len(s) for s in sequences):
    print("Sequence lengths are different!")

# Get all bases.
bases = set(c for seq in sequences for c in seq)

k = 2
kmers = [''.join(p) for p in itertools.product(bases, repeat=k)]
pos_to_freq = []

for i in range(n - k + 1):
    kmers_dict = dict.fromkeys(kmers, 0)
    for s in sequences:
        kmers_dict[s[i:i+k]] += 1/len(sequences)
    pos_to_freq.append(kmers_dict)

for kmer in kmers:
    print("\t".join([kmer] + [str(round(d[kmer], 3)) for d in pos_to_freq]))

```

Output for your example is below (I apologize for the formatting):

```
AA  0   0.143   0   0   0   0.143   0   0
AT  0   0   0   0   0   0.286   0   0
AC  0.143   0   0   0   0   0   0.143   0
AG  0.143   0.143   0.143   0   0   0.143   0.286   0
TA  0   0   0   0   0.571   0   0   0
TT  0.143   0   0   0   0.143   0.143   0   0
TC  0   0   0   0   0.143   0   0   0
TG  0   0.286   0.143   0   0.143   0   0.429   0
CA  0.286   0   0   0   0   0.143   0   0
CT  0.143   0.143   0   0   0   0   0   0
CC  0.143   0   0   0   0   0   0   0.286
CG  0   0.143   0   0   0   0   0   0
GA  0   0   0   0   0   0.143   0   0.429
GT  0   0   0   1.0 0   0   0   0
GC  0   0   0   0   0   0   0.143   0.143
GG  0   0.143   0.714   0   0   0   0   0.143

```

---
Tags: python, sequence-analysis, k-mer
---


# Question
Title: Getting all promoters of a specific gene
I'd like to perform a motif analysis of all promoter regions of the gene ENSMUSG00000020538 in mus musculus. To do so, I wanted to use Biomart :

```
library('biomaRt')
mart <- useMart("ensembl")
mart <- useDataset("mmusculus_gene_ensembl", mart = mart)
all.genes <- getSequence(id = 'ENSMUSG00000020538',
                         type="ensembl_gene_id",
                         seqType="coding_gene_flank",
                         upstream=100,
                         mart=mart)

```

But this provide me only the flanking region of the sepcified gene, not all its promoter regions. How could I do ?

# Answer
*Answer from @ATpoint, converted from comment:*

You need to extract the TSS of all transcripts of that gene, not the gene borders. Query biomaRt for that gene, then set options to extract all transcripts, then get the sequence for the TSS coordinates of those (with a window).
> 0 votes

---
Tags: ensembl, biomart
---


# Question
Title: Splice variant analysis with ensembl genome/annotation/rna. Which files do I use?
I am running splice variant analysis. I wanted to use NCBI genome but the program works better with ensembl. I am a bit confused on primary vs. top level to use as my reference genome. I also do not know which annotation file to use. And I am not seeing a reference rna download. If anyone can help that would be immensely appreciated.

# Answer
> 1 votes
Ensembl provides gene annotation files in both GTF and GFF3 format via the FTP site: http://ftp.ensembl.org/pub/

The GTF and GFF3 files are named according to whether they contain data relating to the primary assembly, or on all top-level sequences including unplaced scaffolds, patches and haplotypes. More information can be found in the README files within each folder. Which genome sequence file and annotation file you use will depend on the details of your study and have been discussed in detail on other posts.

Ensembl also provides reference RNAseq data aligned to the reference genome sequence as BAM files with corresponding .bai index files on the FTP site: http://ftp.ensembl.org/pub/current\_bamcov/homo\_sapiens/genebuild/

More information can be found in the README: http://ftp.ensembl.org/pub/current\_bamcov/homo\_sapiens/genebuild/README

---
Tags: sequence-annotation, transcriptome, reference-genome, ensembl
---


# Question
Title: Find corresponding symbol for gene used in Seurat
I have a Seurat object

```
> pbmc
An object of class Seurat
48608 features across 11709 samples within 2 assays
Active assay: SCT (22220 features, 3000 variable features)
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap
>

```

I want to plot TCRVδ2

```
> FeaturePlot(pbmc,"TCRVδ2")
Error: None of the requested features were found: TCRVδ2 in slot data
In addition: Warning message:
In FetchData(object = object, vars = c(dims, "ident", features),  :
  The following requested variables were not found: TCRVδ2
>

```

This data base is about human TCR genes

https://www.genenames.org/data/genegroup/#!/group/370

But I even can not find this gene here

Do you know what is the corresponding symbol for this gene used in Seurat?

# Answer
`TCRVδ2` is not a gene, but it defines a specific subset of gamma-delta TCRs. You can look at `TRDC` (delta chain) and `TRGC1` and `TRGC2` (gamma chains). `TCRVδ2` should have the combo `TRDC-TRGC1`, while `TCRVδ1` should express `TRDC-TRGC2`. You may need to increase the resolution of your clustering, to avoid clustering the gamma-delta T cells with the NK cells.
> 1 votes

---
Tags: scrnaseq, gene, seurat
---


# Question
Title: Keeping rows with a cut-off in all columns
I have a data frame of raw read counts, genes in rows and samples in columns

```
> dput(head(oc))
structure(list(sample1 = c(13L, 0L, 1L, 42L, 15L, 1137L), sample2 = c(21L,
2L, 20L, 51L, 5L, 203L), sample3 = c(25L, 2L, 6L, 71L, 36L, 1346L
), sample4 = c(59L, 3L, 8L, 58L, 42L, 3428L), sample5 = c(19L,
4L, 8L, 52L, 1L, 201L)), row.names = c("WASH7P", "MIR6859-1",
"CICP27", "WASH9P", "MTND1P23", "MTND2P28"), class = "data.frame")
>

> head(oc)
          sample1 sample2 sample3 sample4 sample5
WASH7P         13      21      25      59      19
MIR6859-1       0       2       2       3       4
CICP27          1      20       6       8       8
WASH9P         42      51      71      58      52
MTND1P23       15       5      36      42       1
MTND2P28     1137     203    1346    3428     201
>

```

I just want to keep genes with at least 10 reads in each sample, I mean even if a gene has less than 10 reads in one sample and the rest of four samples have 10 reads, this gene should be removed. For instance here the second, third and fifth rows should be removed

I tied this but all I obtained zero

```
> df=t(oc)
> cutoff <- 10
>
> # encode the data based on the cutoff
> df[df >= cutoff] <- 1
> df[df < cutoff] <- 0
> View(df)
> rowSums(df)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
>

> apply(df >= 5, 1, sum)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
> apply(df >= 1, 1, sum)
sample1 sample2 sample3 sample4 sample5
      0       0       0       0       0
>

```

Any suggestion please?

# Answer
> 1 votes
`apply(my_df >= 10, 1, sum)` would give you that how many of the samples satisfy your criteria. Basically you `apply()` the function `sum()` to the rows (denoted as `1`). The `sum()` does the trick that @Devon Ryan suggested above.

---
Tags: r, rna-seq
---


# Question
Title: Avoiding the warning 'The following features were omitted as they were not found in the scale.data slot for the SCT assay' in Seurat
I have a Seurat object

I want to plot a list of genes but I got a warning saying most of genes removed

I have tried with both ***my genes of interest*** and ***all variable genes*** but the results is very similar

Does this mean these genes are basically are not expressed in my data?

```
   > t$V1
  [1] "IGH"          "IGHV1-2"      "IGHV1-3"      "IGHV1-12"     "IGHV1-8"
  [6] "IGHV1-14"     "IGHV1-17"     "IGHV1-18"     "IGHV1-24"     "IGHV1-45"
 [11] "IGHV1-46"     "IGHV1-58"     "IGHV1-67"     "IGHV1-68"     "IGHV1-69"
 [16] "IGHV1-69D"    "IGHV1-38-4"   "IGHV1-69-2"   "IGHV2-5"      "IGHV2-10"
 [21] "IGHV2-26"     "IGHV2-70"     "IGHV2-70D"    "IGHV3-6"      "IGHV3-7"
 [26] "IGHV3-9"      "IGHV3-11"     "IGHV3-13"     "IGHV3-15"     "IGHV3-16"
 [31] "IGHV3-19"     "IGHV3-20"     "IGHV3-21"     "IGHV3-22"     "IGHV3-23"
 [36] "IGHV3-25"     "IGHV3-29"     "IGHV3-30"     "IGHV3-30-2"   "IGHV3-30-3"
 [41] "IGHV3-30-5"   "IGHV3-32"     "IGHV3-33"     "IGHV3-33-2"   "IGHV3-35"
 [46] "IGHV3-36"     "IGHV3-37"     "IGHV3-38"     "IGHV3-41"     "IGHV3-42"
 [51] "IGHV3-43"     "IGHV3-43D"    "IGHV3-47"     "IGHV3-48"     "IGHV3-49"
 [56] "IGHV3-50"     "IGHV3-52"     "IGHV3-53"     "IGHV3-54"     "IGHV3-57"
 [61] "IGHV3-60"     "IGHV3-62"     "IGHV3-63"     "IGHV3-64"     "IGHV3-65"
 [66] "IGHV3-66"     "IGHV3-71"     "IGHV3-72"     "IGHV3-73"     "IGHV3-74"
 [71] "IGHV3-75"     "IGHV3-76"     "IGHV3-79"     "IGHV3-38-3"   "IGHV3-69-1"
 [76] "IGHV4-4"      "IGHV4-28"     "IGHV4-30-1"   "IGHV4-30-2"   "IGHV4-30-4"
 [81] "IGHV4-31"     "IGHV4-34"     "IGHV4-39"     "IGHV4-55"     "IGHV4-59"
 [86] "IGHV4-61"     "IGHV4-80"     "IGHV4-38-2"   "IGHV5-51"     "IGHV5-78"
 [91] "IGHV5-10-1"   "IGHV6-1"      "IGHV7-4-1"    "IGHV7-40"     "IGHV7-27"
 [96] "IGHV7-56"     "IGHV7-34-1"   "IGHV7-81"     "IGHVII-1-1"   "IGHVII-15-1"
[101] "IGHVII-20-1"  "IGHVII-22-1"  "IGHVII-26-2"  "IGHVII-28-1"  "IGHVII-30-1"
[106] "IGHVII-30-21" "IGHVII-31-1"  "IGHVII-33-1"  "IGHVII-40-1"  "IGHVII-43-1"
[111] "IGHVII-44-2"  "IGHVII-46-1"  "IGHVII-49-1"  "IGHVII-51-2"  "IGHVII-53-1"
[116] "IGHVII-60-1"  "IGHVII-62-1"  "IGHVII-65-1"  "IGHVII-67-1"  "IGHVII-74-1"
[121] "IGHVII-78-1"  "IGHVIII-2-1"  "IGHVIII-5-1"  "IGHVIII-5-2"  "IGHVIII-11-1"
[126] "IGHVIII-13-1" "IGHVIII-16-1" "IGHVIII-22-2" "IGHVIII-25-1" "IGHVIII-26-1"
[131] "IGHVIII-38-1" "IGHVIII-44"   "IGHVIII-47-1" "IGHV8-51-1"   "IGHVIII-67-2"
[136] "IGHVIII-67-3" "IGHVIII-67-4" "IGHVIII-76-1" "IGHVIII-82"   "IGHVIV-44-1"
[141] "IGHD1-1"      "IGHD1-7"      "IGHD1-14"     "IGHD1-20"     "IGHD1-26"
[146] "IGHD2-2"      "IGHD2-8"      "IGHD2-15"     "IGHD2-21"     "IGHD3-3"
[151] "IGHD3-9"      "IGHD3-10"     "IGHD3-16"     "IGHD3-22"     "IGHD4-4"
[156] "IGHD4-11"     "IGHD4-17"     "IGHD4-23"     "IGHD5-5"      "IGHD5-12"
[161] "IGHD5-18"     "IGHD5-24"     "IGHD6-6"      "IGHD6-13"     "IGHD6-19"
[166] "IGHD6-25"     "IGHD7-27"     "IGHJ1"        "IGHJ1P"       "IGHJ2"
[171] "IGHJ2P"       "IGHJ3"        "IGHJ3P"       "IGHJ4"        "IGHJ5"
[176] "IGHJ6"        "IGHA1"        "IGHA2"        "IGHG1"        "IGHG2"
[181] "IGHG3"        "IGHG4"        "IGHGP"        "IGHD"         "IGHE"
[186] "IGHEP1"       "IGHM"
> pbmc=ScaleData(
+     pbmc,
+     features = t$V1)
Centering and scaling data matrix
  |================================================================================| 100%
> DoHeatmap(pbmc,t$V1,group.by = "Condition")
Warning message:
In DoHeatmap(pbmc, t$V1, group.by = "Condition") :
  The following features were omitted as they were not found in the scale.data slot for the SCT assay: IGHEP1, IGHE, IGHJ6, IGHJ5, IGHJ4, IGHJ3P, IGHJ3, IGHJ2P, IGHJ2, IGHJ1P, IGHJ1, IGHD7-27, IGHD6-25, IGHD6-19, IGHD6-13, IGHD6-6, IGHD5-24, IGHD5-18, IGHD5-12, IGHD5-5, IGHD4-23, IGHD4-17, IGHD4-11, IGHD4-4, IGHD3-22, IGHD3-16, IGHD3-10, IGHD3-9, IGHD3-3, IGHD2-21, IGHD2-15, IGHD2-8, IGHD2-2, IGHD1-26, IGHD1-20, IGHD1-14, IGHD1-7, IGHD1-1, IGHVIV-44-1, IGHVIII-82, IGHVIII-76-1, IGHVIII-67-4, IGHVIII-67-3, IGHVIII-67-2, IGHV8-51-1, IGHVIII-47-1, IGHVIII-44, IGHVIII-38-1, IGHVIII-26-1, IGHVIII-25-1, IGHVIII-22-2, IGHVIII-16-1, IGHVIII-13-1, IGHVIII-11-1, IGHVIII-5-2, IGHVIII-5-1, IGHVIII-2-1, IGHVII-78-1, IGHVII-74-1, IGHVII-67-1, IGHVII-65-1, IGHVII-62-1, IGHVII-60-1, IGHVII-53-1, IGHVII-51-2, IGHVII-49-1, IGHVII-46-1, IGHVII-44-2, IGHVII-43-1, IGHVII-40-1, IGHVII-33-1, IGHVII-31-1, IGHVII-30-21, IGHVII-30-1, IGHVII-28-1, IGHVII-26-2, IGHVII-22-1, IGHVII-20-1, IGHVII-15-1, IGHVII-1-1, IGH [... truncated]
> pbmc=ScaleData(
+     pbmc,
+     features =rownames(pbmc))
Centering and scaling data matrix
  |================================================================================| 100%
> rm(pbmc)
> pbmc=readRDS("pbmc.rds")
> pbmc=ScaleData(
+     pbmc,
+     features =rownames(pbmc))
Centering and scaling data matrix
  |================================================================================| 100%
> DoHeatmap(pbmc,t$V1,group.by = "Condition")
Warning message:
In DoHeatmap(pbmc, t$V1, group.by = "Condition") :
  The following features were omitted as they were not found in the scale.data slot for the SCT assay: IGHEP1, IGHE, IGHJ6, IGHJ5, IGHJ4, IGHJ3P, IGHJ3, IGHJ2P, IGHJ2, IGHJ1P, IGHJ1, IGHD7-27, IGHD6-25, IGHD6-19, IGHD6-13, IGHD6-6, IGHD5-24, IGHD5-18, IGHD5-12, IGHD5-5, IGHD4-23, IGHD4-17, IGHD4-11, IGHD4-4, IGHD3-22, IGHD3-16, IGHD3-10, IGHD3-9, IGHD3-3, IGHD2-21, IGHD2-15, IGHD2-8, IGHD2-2, IGHD1-26, IGHD1-20, IGHD1-14, IGHD1-7, IGHD1-1, IGHVIV-44-1, IGHVIII-82, IGHVIII-76-1, IGHVIII-67-4, IGHVIII-67-3, IGHVIII-67-2, IGHV8-51-1, IGHVIII-47-1, IGHVIII-44, IGHVIII-38-1, IGHVIII-26-1, IGHVIII-25-1, IGHVIII-22-2, IGHVIII-16-1, IGHVIII-13-1, IGHVIII-11-1, IGHVIII-5-2, IGHVIII-5-1, IGHVIII-2-1, IGHVII-78-1, IGHVII-74-1, IGHVII-67-1, IGHVII-65-1, IGHVII-62-1, IGHVII-60-1, IGHVII-53-1, IGHVII-51-2, IGHVII-49-1, IGHVII-46-1, IGHVII-44-2, IGHVII-43-1, IGHVII-40-1, IGHVII-33-1, IGHVII-31-1, IGHVII-30-21, IGHVII-30-1, IGHVII-28-1, IGHVII-26-2, IGHVII-22-1, IGHVII-20-1, IGHVII-15-1, IGHVII-1-1, IGH [... truncated]
>

```

# Answer
> 1 votes
Specify your genes of interest using the `features` parameter of the `ScaleData()` function. If this parameter is not specified, only "variable genes" are scaled and that makes sense, genes that do not demonstrate variability across conditions is not interesting in general. Anyhow, since the `DoHeatmap()` would be looking at the "scaled data slot", it will not be able to find your genes of interest unless they show variation.

---
Tags: scrnaseq, seurat, 10x-genomics, heatmap
---


title: Barplot without alphabetizing x-axis

text:
# Question
Title: Barplot without alphabetizing x-axis
I'm inputting a .csv file and then plotting a barplot based on the colnames but when I go plot the order the names are being alphabetized rather than get the order that is in the .csv file. How to I prevent that from happening?

```
p2 = ggplot(data=clusters_per_sample, aes(x=sample, y=Perc, fill=cluster)) +
  geom_bar(stat="identity")+
  scale_fill_manual(values = getPalette(22))+
  theme_minimal() + labs(y="percentage of total # cells in the sample") + theme(axis.text.x = element_text(angle = 90, face = "bold")) + theme(axis.text.x = element_text(size=12, color="black"))

```

I tried this to get rid of the alphabetizing but my bars are very thin when doing so:

```
p2 = ggplot(data=clusters_per_sample_2, aes(x=sample, y=Perc, fill=cluster)) +
geom_bar(stat="identity")+ scale_x_discrete(limits=clusters_per_sample_2$sample) +
heme_minimal() + labs(y="percentage of total # cells in the sample") + theme(axis.text.x = element_text(angle = 90, face = "bold")) + theme(axis.text.x = element_text(size=12, color="black"))

```

Any suggestions?

# Answer
> 0 votes
ggplot2 uses factors on the underlying data frame. To have the data plotted in a custom order, apply a custom ordered factor on the underlying dataset.

For example, if `clusters_per_sample$sample` is your custom-ordered vector, do this:

```
clusters_per_sample$sample <- factor(clusters_per_sample$sample, levels = unique(clusters_per_sample$sample), ordered = TRUE)

```

Now when you plot `sample` as an aesthetic, ggplot has a factor component to work with already, so it won't try and make a new one (this new factor would have been in alphabetical order, which we circumvent with the `ordered=TRUE` parameter in our custom factor).

---
Tags: r, ggplot2, graphs
---


# Question
Title: sctransform - mitochondrial expression filtering / transformation
I have a Seurat object which has a high expression of `mitochondrial genes`

```
> pbmc1
An object of class Seurat
36601 features across 18338 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

```

I want to use `SCTransform in Seurat`, I don't know if still I should define the mitochondrial percentage or not, if not, how SCTransform function knows which percentage suits for my data?

I used SCTransform without clipping `MT genes` and I noticed I have the same number of cells

```
> pbmc_SCTransform <- SCTransform(pbmc, method = "glmGamPoi", vars.to.regress = "percent.mt", verbose = FALSE)
> pbmc_SCTransform
An object of class Seurat
58629 features across 18338 samples within 2 assays
Active assay: SCT (22028 features, 3000 variable features)
 1 other assay present: RNA

```

Is it fine or I should clip `MT genes` before running SCTransform function?

I am not too sure about the input composition but from 10X websummary two of my samples likely suffer from poor sample prep or poor beginning cell health. The gene expression median genes per cell are 41 and 50 and differential expression show only mitochondrial (MT) and MALAT1 genes. This is indicative of subpar lysis conditions.

# Answer
Mitochondrial gene filtering is done on a case-by-case basis, which is why there's a bit of discussion about choosing appropriate thresholds in the Seurat PBMC3k vignette.

In your specific case, there are a lot of cells with very high mitochondrial gene content. It would be a good idea to work out why that is before doing any more filtering.

However, you're asking specifically about `sctransform`, and the SCtransform vignette does account for mitochondrial gene content in the transformation:

> During normalization, we can also remove confounding sources of variation, for example, mitochondrial mapping percentage

```
# store mitochondrial percentage in object meta data
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")

# run sctransform
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)

```

There's no filtering in that vignette, but there's also no mention of how clean the data was initially. It's still a good idea to look at your data for outliers, try to understand why they are outliers, and have a think about whether or not to exclude them.
> 1 votes

---
Tags: scrnaseq, seurat, 10x-genomics
---


title: R limma alternatives in Python

text:
# Question
Title: R limma alternatives in Python
The R package `limma` is ideal to perform differential expression analysis. Is there any `limma` alternative in Python?

I'm trying to use `statsmodels` and `scikitlearn` in conjunction with some other cool tools (such as `pycombat`) to get `limma`-like workflows. However, I'm missing many capabilities from the `limma` toolbox.

# Answer
**August 2021**

To date, it seems that the response to this question is **"No"**. However, there is a GitHub repository called `edgePy`, aiming to *"become an implementation of `edgeR` for differential expression analysis in the Python language"* (`edgeR` is somewhat similar to `limma`).

PS: Many thanks to @ATpoint90 (the Twitter "hive mind" summoner) and @sridnona6 (`edgePy` developer).
> 2 votes

---
Tags: r, python, differential-expression, limma
---


# Question
Title: Comparing between different set of clusters produced from the same network
I have a set of eight gene clusters which I obtained from clustering a network. However, I had to update the genes in the network at a later stage which lead to certain deletions (genes not present in the new gene list). I did the clustering again and this time there were a total of 9 clusters. Is there some way in which I could compare the earlier eight clusters with the new nine ones to effectively say that "clustering has happened along the same lines". It would have been easier to do an eight to eight comparison checking from the largest to the smallest cluster. However, since I have nine clusters now, I am kinda in a fix as how to compare. I used the same clustering algorithm to cluster the network in both the cases

# Answer
*Answer from @llrs, converted from comment:*

Perhaps I'm not understanding correctly what you are doing. But these stability tests are used when the number of clusters is unknown, so you test the resolution, or other metrics for a range of numbers of clusters. Then you choose the best one (either the one that is better by the score, or the one that increasing the number of clusters don't improve your metric much more)
> 0 votes

# Answer
If you have both cluster assignments, you can check for consistency of assignment by looking at pairwise counts using `table` in R:

```
cons.mat <- table(clusters.1, clusters.2)

```

The result of this will be a matrix showing the number of cells / samples that have the observed cluster assignments via each method. If the two cluster assignments are similar, then samples from one cluster via one method should be mostly assigned to one cluster in the other method. You can check this by comparing table cell values with their row/column sum:

```
## show consistency along the first dimension
apply(cons.mat, 1, function(x){max(x / sum(x))})
## show consistency along the second dimension
apply(cons.mat, 2, function(x){max(x / sum(x))})

```
> -1 votes

---
Tags: networks, clustering
---


# Question
Title: Hisat2 compatibility for long reads (Pacbio)
I am working with a high performance cluster computer containing 112 threads. I am trying to align PacBio transcriptome reads against the genome to count the gene number. For pair end read i used the following workflow:

```
#convert gff to gtf
/home/software/cufflinks-2.2.1/gffread xxx.gff -T -o xxx.gtf
#build index
/home/software/hisat2-2.2.1/hisat2_extract_exons.py xxx.gtf > xxx.exon
/home/software/hisat2-2.2.1/hisat2_extract_splice_sites.py xxx.gtf > xxx.ss
/home/software/hisat2-2.2.1/hisat2-build -p 20 xxx.fa --ss xxx.ss --exon xxx.exon xxx_tran
#align reads
/home/software/hisat2-2.2.1/hisat2 -p 80 -x ./xxx_tran -1 R1.fq -2 R2.fq -S gonad.sam

#convert to bam and sort by read name
samtools view -bS gonad.sam -o gonad.bam

samtools sort -@ 80 gonad.bam -o gonad.sorted.bam

#count gene
htseq-count -f bam --strand=no gonad.sorted.bam xxx.all.gtf > gonad-counts.txt

```

But in case of Long read i am not sure what i need to do during the alignment stage. I have tried with minimap2 like as following:

```
#align reads
/home/software/minimap2-2.17_x64-linux/minimap2 -ax splice:hq -uf yyy.fa gill.fq > gill.sam

```

in following, i stucked in counting stage:

```
#count gene
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt
error txt:
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

```

My question is that during counting it requires sorted bam and genome gtf file that was generated succesfully during previous staps, why it is getting failed even though the same pipeline is absolutely fine for pair end reads? Hisat2 can align long read? Is there any other way do that? As a novice i will be always grateful to you. Thanks.

Update: I tried generating an index:

```
# generate index
samtools index gill.sorted.bam

```

And the output was a .bai file like "gill.sorted.bam.bai". And I tried again to count as written above. But then I found the following error:

```
100000 GFF lines processed.
200000 GFF lines processed.
300000 GFF lines processed.
400000 GFF lines processed.
438395 GFF lines processed.
Error occured when processing input (record #79 in file gill.sorted.bam):
  'NoneType' object has no attribute 'encode'
  [Exception type: AttributeError, raised in _HTSeq.pyx:1379]

```

Is there any way?

Update 2:

Here is what the BAM file looks like:

```
SRR6418045.89565    2048    PYXC01000001.1  9893    60  5794H11M1I42M1I3M2I19M1D20M1I15M1D20M3D43M1I11M3D8M1I19M3I24M1I7M1I4M1I5M1D6M1D14M1I45M1D2M1D3M1D4M1I8M1I25M1D21M1D62M1D3M1D25M2D21M2I30M1I25M1D2M1D5M2I4M1D13M1D12M1D10M1I59M21564H    *   0   0   GAATACATAAACTGGATAGCGCACTGTTAATATTTACAGTGAAGAGAGTCTCTTCGAATTTTCCTCCTCAATCATCTTATTAGAAAGATGTTGAGTGGGCAGCAGCCAAATAAAATACTGTTCAGTGATCACTTCCAATGATGCTGTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGCCTAATTTAAAAAACTAATGTAAATTATTTAAAAGTTTTTTTTTTTTTTATGAACTTTTTAGCCACATTGAATGGCTTCAGAATATCTACAGCCTAATGTATTTTCTTAGTAGCCTAATGAAGAAAGAATATTACCTTAAGTTACCATAGATTGACGCCAAAACTCGTCAAATCTCAAACCCCTAACATTTGTTCAACTGTTGTATACATGCTCATACACAGTGTTTCTTGTGTTCCTTAACCCACCGGGTAAAAAGTGGAGGTCCATGTCTGATACGTAGTTATCAATGTACTGTTGCCCATTTAGAGAGCACTTTGGTTTTCCCTCACCAGATAAGGTTTTTTCTCACTTTCAGTGTAAACCTGAGACTTGTTTGCACTTAATATTGGGGATTGTGTAACTCTTATTCAAGGAATTGCCTGTTTTCCAATTTGGTCTGATTACACTAACCTTTGTGTAAAATTATGGTTCAAATAAACTGTATCACA    -)/////,+$/"/..+//-/-+/.-/-.*,./.,)//'//.+/./..-/./../"/*/.),#)*.)-+)*-(///&.//.*'/,*////*-,./+.++#",.///,/,*,/-,)/,/)/.-)/.///./-..-*++)-(/*,../..)+//.////+*...////./...(././,/.,*.-'*/)(.-*)&(#'&,(#'....').,,$*')+#&,,)%%%$&,)*)*(*,..,'/-.+*)*,)+',-).-)*//,-",-+./)///'-.//,-+,.*//.-/%,**)-..../.*/,*'/.*//,*/.,(&.,+-)-'-*/*-/.-/-&+$&,'.'.&,.*,)/-."/,.,*/-//.,)/*,%.+*)/+.,--,./*--/.-.//'///)/--/$&//-//.//.-)/),,./-)+&../*,*,..-/+,)/,,,-*,/-,/.(.(///./(////./&-./-/+//*...////,&%+$&%#%$$--.././/.-+',-,#(.&%(-/,*/.///*..--++*)/-.*..*)*(,$*/--*.+./'//...-..)-.&.-)&))+.-,&#'-'--..'/'-//',//-*.)/+-*....%../.'(/*-$./,(.(////.#/..+/+%/'-'%...//*$+#--+/&-,-/,*-.-*/-.&../-,%/    NM:i:55 ms:i:343    AS:i:343    nn:i:0  ts:A:+  tp:A:P  cm:i:77 s1:i:460    s2:i:0  de:f:0.0656 SA:Z:PYXC01000001.1,9893,+,4147S672M2I23209S,60,37;PYXC01000001.1,9893,+,847S669M26I26488S,60,42;PYXC01000001.1,9893,+,26975S673M3I379S,60,44;PYXC01000001.1,9910,-,10980S654M9I16387S,60,41;PYXC01000001.1,9900,-,20736S661M9I6624S,60,43;PYXC01000001.1,9893,+,7409S673M15I19933S,60,49;PYXC01000001.1,9893,-,27284S667M12I67S,60,52;PYXC01000001.1,9893,-,19099S673M3I8255S,60,48;PYXC01000001.1,9893,-,7656S673M1652I18049S,4,1688;PYXC01000001.1,9908,+,2517S649M17I24847S,60,51;PYXC01000001.1,9893,+,23715S672M1D3643S,60,51;PYXC01000001.1,9893,-,25633S669M22I1706S,60,53;PYXC01000001.1,9893,-,17479S668M1I9882S,60,53;PYXC01000001.1,9893,-,14214S673M25I13118S,22,55;PYXC01000001.1,9899,+,13933S663M13I13421S,60,55;PYXC01000001.1,9900,-,22352S658M21I4999S,60,55;PYXC01000001.1,9901,+,15556S665M7I11802S,60,60;PYXC01000001.1,9893,+,12284S668M9I15069S,60,58;PYXC01000001.1,9909,-,24008S657M16I3349S,60,57;PYXC01000001.1,9902,+,25358S661M10I2001S,60,60;PYXC01000001.1,9893,+,10659S671M2D16700S,60,60;PYXC01000001.1,9893,+,9039S668M10I18313S,60,51;PYXC01000001.1,9893,-,2790S673M21I24546S,60,65;PYXC01000001.1,9899,+,22114S669M8I5239S,60,59;PYXC01000001.1,9893,-,6025S673M17I21315S,60,64;PYXC01000001.1,9908,-,12604S645M13I14768S,60,64;PYXC01000001.1,9893,-,4418S673M3I22936S,60,71;PYXC01000001.1,9893,-,1162S673M15I26180S,60,70;PYXC01000001.1,9898,+,17178S641M1I10210S,60,63;PYXC01000001.1,9893,+,20477S673M13I6867S,60,67;PYXC01000001.1,10390,-,64S160M2D27806S,60,8;   rl:i:2011
SRR6418045.89565    2048    PYXC01000001.1  9893    60  12284H12M1D27M1D9M1I15M1I26M1D19M2I16M2I2M1I17M1I32M1I31M2D8M3I25M1I6M1I19M2I12M1D32M2D15M1I6M2I7M1D4M2D17M1I10M1I32M1I36M1D4M1D2M1D20M1I59M1I11M1D10M1D10M1I24M1I9M2D2M1D28M1I24M1I11M15069H   *   0   0   GAATACATAAATGATAGCGCACTGTTAATATTTACAGTGAGAGAGTCTACTTGAATTCCTCCTCAAATCATCTTATTTAGAAAGATGTTGATGGGAGCAACCAAATAAAAATCTACTGTTCAGTGATCCAACTTTCACAAAATGATGCTGTTGCACTTACAGAAGCATGTATGGTACACAGTGTTATGACTATTTTAATAAAAAACAATGTAAATTTTAAAAGGTTTTTTTTTTTTTATGAACTTTTTAGCAACATTGAATGGTTCAGAAATATCTTAGCAAGCCTAATGTATTCTTAGTAGCCTAATGAAGAAAGAATATTACTAAGTTACCAGTAGGAATCTGAACGCCAAAATGTCATCTCAAACCCCTAACATTTTGGTTCAATCTGTTGTATACATGCATCATACACAGTGTTTCTTTGTGTTCCTTTAACCACCGGGTAAAAAGTGGAGGTCATGCTATTACTGTAGTTATCAATGTAACTGTTGCCCACTTTTAGAGAGCACTTTGGTTTTCTCACCAGATAAGGTTTTTTCTCACTTTTTATGTAAACTGAGACTTGTTGCACTGTACTATATTGAATATGTGTAACTCTTATTTTCAAGGATGCCTGTTTTCCATTTGGTCTATTAACACTAACACTTTGTGTAAACTTATGGTTCAAAATAAACTGT   -*//..//-*./%/,/./,,///.,/+//*/,////..,*..+/..-,(*.,.%&".*./%-..((//$///...)+-&.+'/////.+*$(),,../$.'.,)//-++$.."//,-.,/////'//,/$...*//./-()./+//,///,...%/-//,*.*/.//&/.//...///./..(-.///#/.+$(#.'..((%*(.)(-+-..%.*,((,,(..(-&+*)*******..+/$-..,++..--$-..'-.)//..-)/.,,*...//./.#.&/)*/)*(/./*,.&-(+////*/.).).+//-*//*)/*.*..-(//-/&,/-//..)*//.+.$/.',-+,(+(&//..//.,).%')./*///.--)--*.'(/"//...//.././/-)'(/,'..)%/./..&..$*,+-.-/*/++-%.+--+.),-,#,&(//.-,).-/./-.'./..///....-(*-/)/.#-$.//.%$*(*..())+/--/.../..(,+,..*).-/#.*'.+/-(/...&'+)/--.+'-))&$,'//*'/.-/(/.).,/,-,*.,/.*.".)..)./*/*/+/.'/%/&-),/..,)).(/--/&,*/..,*)-#+/.,-,./-**./%.//+/*)().(+-,,-.%'//../--.--/)*#..,*///.-   NM:i:58 ms:i:325    AS:i:325    nn:i:0  ts:A:+  tp:A:P  cm:i:65 s1:i:397    s

```

# Answer
This error:

```
[E::idx_find_and_load] Could not retrieve index file for 'gill.sorted.bam'

```

Means that you have not generated an index for a BAM file. This needs to be done for most downstream processing of BAM files:

```
# generate index
samtools index gill.sorted.bam
# count genes
htseq-count -f bam --strand=no gill.sorted.bam yyy.gtf > gill-counts.txt

```

Update (after second error):

I've got no idea what's going wrong; this is not an error that I would expect to see, and doesn't look related to PacBio input. I wonder if the error is related to the number of threads used in the sort command, given that it happens at record 79, and you've specified 80 threads. The next things that I would do would be to look at the input files around where the error is to try to work out what's going wrong:

```
samtools view gill.sorted.bam | head -n 85 | tail -n 10

```
> 1 votes

---
Tags: gene, transcriptome
---


title: Slow bwa shm "preloading" into memory

text:
# Question
Title: Slow bwa shm "preloading" into memory
Coming from devops background, working on automating usual genomics pipelines, I am benchmarking existing pipelines in order to choose cloud resources properly.

Past week most of work was on testing `bwa mem`.

I've noticed that `bwa` takes about 10 minutes to start processing data. When monitoring `bwa` resource usage I've noticed that RAM is slowly building up to about 6GB. That is amount of reference indexed data.

After bwa loads the reference first time, any subsequent run is fast. But why this preloading is happening so slow? My reference data set is already indexed with `bwa index`.

When I monitor disk activity with `iotop` I noticed that `disk read` is only `5MB/s`, when I run `bwa mem` alignment (only first time), which is very slow. Performance of the drive is lot more than that and a simple `copy` test confirms that.

Is it a "limitation" of `bwa`, is this slow preloading expected to be like that and I should just accept it and move on forward or there is something I don't understand and maybe someone can help me understand what could I do?

The same "slow preload" happens when I run `bwa shm /path/to/ref`. Then `bwa mem` starts immediately. But again slow preload is something that is bothering me and I am trying to understand what can I do about it.

# Answer
> 1 votes
I'll add an answer to my own question, which was asked out of lack of experience with cloud (AWS) EC2 service.

The way I was testing `bwa` is from a custom AMI image created for our pipelines.

The image has "preloaded" reference data that we are using for alignement so that we spare some time downloading reference data from `S3` and spare some `S3` api calls since we have a lot of files used as reference data.

Simply our image is single 200GB EBS volume used with `AWS ECS AMI` where we stored our reference data and created a custom image for Batch out of it. Then we use that AMI for AWS Batch jobs.

Once I used `iotop` (after few days of testing) to monitor disk activity I've noticed that `bwa` loads index into memory very slow, cca `5MB/s`. Syncing files from S3 was faster then using this `snapshot` approach.

For benchmarking I've used the same AMI manually launched. In the AWS documentation there is part which explains why and how you should optimize EBS when AMI is launched from a snapshot. https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ebs-initialize.html

After spining up fresh EC2, syncing reference data and running `bwa` everything was clear. My EBS volume was slow and had to be initialized first.

That is all. My lack of EBS storage knowledge was the issue.

This helped me understand that I should not use custom AMI with ref data stored and will probably use EFS or simply syncing the data from S3 inside Docker entrypoint.

# Answer
> 0 votes
*Answer from @user438383, converted from comment:*

BWA has slow indexing due to the indexing size; this can be improved by using bwa-mem2:

https://github.com/bwa-mem2/bwa-mem2

> We are happy to announce that the index size on disk is down by 8 times and in memory by 4 times due to moving to only one type of FM-index (2bit.64 instead of 2bit.64 and 8bit.32) and 8x compression of suffix array. For example, for human genome, index size on disk is down to ~10GB from ~80GB and memory footprint is down to ~10GB from ~40GB. ***There is a substantial reduction in index IO time due to the reduction*** and hardly any performance impact on read mapping.

---
Tags: bwa
---


# Question
Title: calculating phylogenetic tree from pairwise distances
I have a file with a 3-column list of pairwise distances:

```
A B distance
A C distance
B C distance
...

```

And I would like to calculate a phylogenetic tree for all the nodes in the list. I've seen `phylip neighbor` and similar options, but it seems like the length of the ids is a limit, which I rather not have to deal with. What would be a straightforward way of getting a newick tree output out of an input as above?

# Answer
> 3 votes
I believe there are a number of ways to construct a tree metric from a distance metric. A very straightforward method, neighbor joining, is available in the Sciki-Bio package. A less straightforward option with more freedom is by using the scipy.cluster.hierarchy module to obtain a linkage matrix, then using the `to_tree()` method to obtain a scipy `Tree` object, and finally writing a script to convert a scipy Tree object to a newick string.

I've copied the example from scikit-bio below:

```
>>> from skbio import DistanceMatrix
>>> from skbio.tree import nj
>>>
>>> data = [[0,  5,  9,  9,  8],
...         [5,  0, 10, 10,  9],
...         [9, 10,  0,  8,  7],
...         [9, 10,  8,  0,  3],
...         [8,  9,  7,  3,  0]]
>>> ids = list('abcde')
>>>
>>> dm = DistanceMatrix(data, ids)
>>> # Or from a tsv file with header, no column 0
>>> reader = csv.reader(open(args.inputfile), delimiter="\t")
>>> x = list(reader)
>>> dm = DistanceMatrix(x[1:],x[0])
>>>
>>> newick_str = nj(dm, result_constructor=str)
>>> print(newick_str[:55], "...")
(d:2.000000, (c:4.000000, (b:3.000000, a:2.000000):3.00 ...

```
* Note: Both methods require you to supply the input as a numpy distance matrix

---
Tags: phylogeny
---


# Question
Title: Using grep to get lines from a file (.tsv) that contain specific samples listed in a sample file (.tsv)
So I have a file, let's call it pcawg.tsv. It's formatted like the one down below.

```
Tumor sample ID  Normal sample ID  Mechanism
a                1                 deletion
b                2                 ecDNA
c                3                 ecDNA

```

I have a file of tumor sample ID's with the ecDNA mechanism called ecDNA\_samples.tsv. It's formatted like this:

```
Tumor sample ID
a
b

```

To filter and extract the lines of `pcawg.tsv` that contain the tumor sample ID's in `ecDNA_samples.tsv` into a new file, I used this line of code below.

```
grep -f ecDNA_samples.tsv pcawg.tsv > ecDNA.vcf.tsv

```

Would this be the right code to use?

# Answer
This may help:

```
cat ecDNA_samples.tsv | xargs -I % grep '^%\s' pcawg.tsv

```

The command above executes `grep '^<line>\s' pcawg.tsv`, for each line of ecDNA\_samples.tsv. Which matches the lines taken from ecDNA\_samples.tsv to the sequences of characters between line start and a white-space in pcawg.tsv.
> 1 votes

# Answer
Just in case an R solution is acceptable as well, this should be fairly robust code because it properly parses the TSV files and matches up column names:

```
#!/usr/bin/env Rscript

library(tidyverse)
read_tsv("pcawg.tsv") %>%
   inner_join(read_tsv("ecDNA_samples.tsv")) %>%
   write_tsv("filtered_pcawg.tsv")

```
> 0 votes

---
Tags: linux, cancer, filtering, awk, grep
---


# Question
Title: How to read and edit a FASTA file with python using regular expression?
I have a file containing protein sequence information of 100 odd proteins. Each sequence starts with a header that looks like this:

```
>sp|Q9UDW1|QCR9_HUMAN Cytochrome b-c1 complex subunit 9 OS=Homo sapiens OX=9606 GN=UQCR10 PE=1 SV=3

```

I want to edit these header to look like this:

```
>Q9UDW1

```

I have tried the following code:

```
import re
with open ('WT.fasta', 'r+') as f:
    contents = f.readlines()
    for i in range(len(contents)-1):
        target = re.compile(r'>..\|(.*)\|.*', contents[i])
        target.sub(r'\1', contents[i])

```

I think I am not being able to call the right methods to do this. I have been trying to read the documentation: https://docs.python.org/3/library/re.html. I am not been able to understand it and feeling overwhelmed by the details. I would really appreciate if anyone could help or point me in the right direction. Thanks!

# Answer
In general, I usually recommend Biopython for any sort of fasta parsing, but the code below should work as well.

```
with open ('WT.fasta', 'r+') as in_f, open("WT_out.fasta", 'w') as out_f:
    for line in in_f:
        if line[0] == ">":
            out_f.write(">" + line.split('|')[1] + "\n")
        else:
            out_f.write(line)

```
> 3 votes

# Answer
You could do this with `awk` by processing the header:

```
$ awk -v FS="|" '{ if ($0 ~ /^>/) { printf(">%s\n", $2) } else print $0 }' in.fa > out.fa

```

For example:

```
% echo -e ">sp|Q9UDW1|QCR9_HUMAN Cytochrome b-c1 complex subunit 9 OS=Homo sapiens OX=9606 GN=UQCR10 PE=1 SV=3\nACCGT" | awk -v FS="|" '{ if ($0 ~ /^>/) { printf(">%s\n", $2) } else print $0 }'
>Q9UDW1
ACCGT

```
> 3 votes

---
Tags: python
---


title: DESEQ2 Question about results()

text:
# Question
Title: DESEQ2 Question about results()
I am currently learning to perform Differential Analysis via DESEQ2 R Package, and I believe I've made progress, able to format the data correctly \[maybe\] for DDS(). When I run the results function to see the output, the data seems fine, as can be seen below:

```
log2 fold change (MLE): Group RA vs C  Wald test p-value: Group RA vs
C  DataFrame with 6 rows and 6 columns
       baseMean log2FoldChange     lfcSE      stat    pvalue      padj
      <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric> Gnai3   95.1978      0.3458366  0.275366  1.255916  0.209146  0.999383
Cdc45   21.9942     -0.0545177  0.346199 -0.157475  0.874870  0.999383
H19     29.0133     -1.4954220  1.516390 -0.986173  0.324048  0.999383
Narf    16.6967      0.2571059  0.370880  0.693232  0.488164  0.999383
Cav2    63.6356      0.5060019  0.354969  1.425480  0.154018  0.999383
Klf6    65.8469      0.2851175  0.265755  1.072858  0.283335  0.999383

```

However, when I run summary() or any other function, they do not seem to behave correctly, not counting any stats. Same with MAplot not coloring genes that should be marked red:

```
out of 13184 with nonzero total read count adjusted p-value < 0.1 LFC
> 0 (up)       : 0, 0% LFC < 0 (down)     : 0, 0% outliers [1]       : 150, 1.1% low counts [2]     : 0, 0% (mean count < 1)

```

I'm confused as to why this is occurring and was wondering if anyone had ideas.

EDIT: Added More Code for context

```
sergey <- (read.csv("cleanmatrix25+.csv", sep = "\t", header = T,
    check.names = F, row.names = "Gene_names"))

meta <-(read.csv("metadata25+.csv", sep = "\t", row.names = 1))

dds <- DESeqDataSetFromMatrix(countData = sergey, colData = meta,
    design = ~Group)

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds$Group <- relevel(dds$Group, ref = "C")

dds <- DESeq(dds)

res <- results(dds)

head(res)

summary(res)

plotMA(res, ylim=c(-2,2))

```

EDIT2:

I sorted result by p-value, so all of these should be red on the following MAplot due to having a p-value \< 0.1

# Answer
You do not have any significantly differential genes after FDR correction according to padj, simple as that, which means that either there are none or your sample size is too small (underpowered) given the observed variability of your data. As the MA-plot shows some logFCs that could qualify as large enough to be called DEG it probably comes down to a too small sample size and quite some variability between the replicates of the groups. A PCA can help diagnosing whether there is simply lots of per-group variability or other factors such as batch effects that need to be taken care of. Please see the DESeq2 vignette for code examples or Bioc packages such as `PCAtools` for this.

Edit: `plotMA` by the way colors based on padj, not p-value alone.
> 2 votes

---
Tags: r, deseq2
---


title: Seurat - single cell analysis

text:
# Question
Title: Seurat - single cell analysis
I have been doing single cell analysis based on seurat and I realized after each clustering and reclustering when I run SCTransform, in my code there is vars.to.regress = percent.mt and I was running this everytime. I removed MT- contaminated cells before.

My clustering seems fine though.

Do you think this is very critical thing and I have to analyze everything from scratch?

# Answer
I think it would be a good idea to at least look at the distribution of mitochondrial gene fraction in your cells. If they mostly have low fractions (e.g. less than 15%), then I wouldn't expect the clustering to be influenced too much by them.

With SCTransform regressing on the percent.mt, I think it's okay to do post-transform filtering, but it'd be a good idea to exclude any cells that have high mitochondrial fraction prior to doing clustering (unless, for whatever reason, those are the cells you're interested in).
> 1 votes

---
Tags: seurat
---


# Question
Title: Denovo, Stacks: Getting an “ambiguous redirect” error
I tried running the following command

```
/public/home/WHC_BBT/Softwares/bin/bin/denovo_map.pl --samples $RawData -O $Popmap -o $fullM2n3 -M 2 -n 3 -m 3 -T 15 &> $logfile

```

This is the error I got

> /opt/gridview//pbs/dispatcher/mom\_priv/jobs/14923.admin.SC: line 6: $logfile: ambiguous redirect

I need help on how to resolve this.

Thanks.

# Answer
This most likely means your variable `$logfile` does not exist.
> 1 votes

# Answer
At a guess, the variable `$logfile` is not properly defined (i.e. it hasn't been defined at all). What is the output of the following command:

```
echo "logfile: ${logfile}"

```
> 0 votes

---
Tags: snp, assembly, radseq
---


title: Adjusted P-value from DESeq2

text:
# Question
Title: Adjusted P-value from DESeq2
What does the padj really mean? I believe it's a BH method, where it attempts to quantify the likelihood of a "false positive" but I am not sure. Also, what does it mean if the resulting calculations nearly all have the same padj values despite having different p-values?

# Answer
> 1 votes
If we look at documentation of DESEQ2 and search "adjusted p-value", we find the section "Multiple test correction". In this section they discuss Benjamini-Hochberg ("BH") false discovery rate (FDR) correction procedures:

> FDR/Benjamini-Hochberg: Benjamini and Hochberg (1995) defined the concept of FDR and created an algorithm to control the expected FDR below a specified level given a list of independent p-values. **An interpretation of the BH method for controlling the FDR is implemented in DESeq2 in which we rank the genes by p-value, then multiply each ranked p-value by m/rank.**

Earlier, they write:

> m = total number of tests

They further add:

> So what does FDR \< 0.05 mean? By setting the FDR cutoff to \< 0.05, we’re saying that the proportion of false positives we expect amongst our differentially expressed genes is 5%. For example, if you call 500 genes as differentially expressed with an FDR cutoff of 0.05, you expect 25 of them to be false positives.

I would suggest reading Benjamini and Hochberg (1995) to learn more details if desired.

---
Tags: statistics, deseq2, p-values
---


# Question
Title: PyMol: Output Polar Contacts Between Chains to Text File
I am new to PyMOL but have a very specific task that I need to do.

I have a PDB structure file of protein homo-oligomer, and I want to use PyMol to determine polar contacts between chain A and all other chains.

**However, I not only want to visualize these polar contacts (which is easy to do via PyMOL GUI), but I also want to be able to store/save pairs of interacting residues** in a dictionary or whatever may be the most convenient data structure to do this. The end goal to create a sort of "interaction map" in a csv or text file where each interacting pair of residues has its own line.

I would ideally like to do this via a Python script where I call the PyMol API or through a PyMol script - whichever is easier.

I have seen a similar question on Stack Overflow but found it unhelpful because I not only want to select and visualize these polar contacts but also want to save/record them in a list or dictionary (and then to a text or csv file), which the above question does not address.

# Answer
So you can install PyMOL as a standalone, but you can also install it as a bona fide Python 3 module via conda:

```
conda install -c schrodinger -c conda-forge pymol-bundle

```

It can be installed in other ways —apt-get or brew or even compiled, but the latter is excruciatingly painful.

In your python notebook you can do:

```
import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('structure_file.pdb', 'foo') # file
    # pymol.cmd.fetch('PDB_CODE') # PDB code
    # pymol.cmd.read_pdbstr(foo_str, 'foo') # string
    ...

```

where the ellipsis can be any API cmd.

This context manager approach is to avoid the problem of namespace pollution as the session (e.g. loading structures without removing the old ones etc.) is the same but you don't see what is happening. Also PyMOL() is parallel the parallel mode (so mutagenesis wizard does not work), while in singleton mode (`pymol2.SingletonPyMOL`), you get issues if you try to re-initialise etc.

Anyway. To "select" all interacting atoms at the interface you can do:

```
n = pymol.cmd.select('chain A around 3')

```

`Expand` includes the selection. `Around` does not.

This adds a selection to the session. So what we actually want is a pythonic instance to play with. So:

atoms = pymol.cmd.get\_model('chain A around 3')

where atoms is a list of `chempy.Atom` objects. With attributes such as resi, resn, chain etc. Do dir on an atom for more...

One attribute is `coord`, which contains a length 3 list of floats.

NB. The above selection gives not a full residue, whereas 'byres chain A around 3' does. So if you want to make a list of all interacting residues:

```
interactors = [{'chain': atom.chain, 'index': atom.resi, 'number': atom.resn} for atom in pymol.cmd.get_model('name CA and byres chain A around 3')]

```

An iteration can give the atom to atom distances.

Do note that this is all proximity, so you will have to screen out hydrogens and carbons to get actual polar contacts.
> 2 votes

---
Tags: python, protein-structure, visualization, pymol, protein-protein-interaction
---


# Question
Title: How to check if a fasta file and a GTF file fit and form a valid pair?
I'd like to know if there is a simple way to check the concordance between a GTF file and a fasta file (that I use as mapping reference) ?

It may be a dumb question, but I suspect discrepancies between my 2 files and I don't want to check this manually

# Answer
> 1 votes
*Answer from @ATpoint, converted from comment:*

Check whether chromosome identifiers are the same and both files contain the same chromosomes. You can also check that the end coordinate of the GTF never exceeds the length of the respective chromosome. Maybe you can also check whether the first three nucleotides of every CDS indeed start with a start codon. If this is all true then I guess I would trust this pair.

# Answer
> 0 votes
You might use `agat_sq_filter_feature_from_fasta.pl` from AGAT

```
 DESCRIPTION
 This script is a kind of annotation filter by sequence name. It goes through the gff annotation features and remove those
 that are not linked to a sequence from the fasta file provided. The
 match between sequence name in the fasta file and the 1st column of
 the gff3 file is case sensitive.

```

---
Tags: fasta, gtf
---


# Question
Title: Limma decideTests function: what kind of multiple hypothesis testing correction does parameter "method" involve?
What kind of multiple hypothesis testing correction does `method="global"` do in Limma's `decideTests` function? According to the documentation:

> method="global" will treat the entire matrix of t-statistics as a single vector of unrelated tests.

and

> method="global" is useful when it is important that the same t-statistic cutoff should correspond to statistical significance for all the contrasts.

What does this really mean? Does it use a Bonferroni correction?

# Answer
More context from the docs, looking at the arguments to `decideTests()`:

> method: character string specifying how genes and contrasts are to be combined in the multiple testing scheme. Choices are "separate", "global", "hierarchical" or "nestedF".

> adjust.method: character string specifying p-value adjustment method. Possible values are "none", "BH", "fdr" (equivalent to "BH"), "BY" and "holm". See p.adjust for details.

These are two different options.

One option (`adjust.method`) chooses the multiple test correction, for which the default is BH/FDR.

The other option (`method`) decides how the multiple test correction is applied, for example can it be applied separately to different genes or across contrasts. `method = global` in this context just means that every single test statistic is considered at the same time. An alternative to this would be to consider a nested design.

More intuitively, `method = separate` applies the correction to each column of your matrix independently, such that in a $m x n$ matrix of statistics, the "number of tests" is $m\*n$ for `global` and $m$ for `separate`.

For understanding the effect of the number of tests on FDR, I suggest consulting the FDR wiki page. Basically, using a `method` other than `global` will increase power.
> 1 votes

---
Tags: rna-seq, limma
---


# Question
Title: Best practice for running GATK VQSR on X chromosome
According to GATK best practice, it is recommended that different VQSR models be built for SNPs and INDELs, because the annotations for high-quality SNPs and INDELs are systematically different (if I understand it correctly). Since annotations for good variants on autosomes could be different from those on X chromosome, e.g., DP for good variants on X chromosome could be substantially smaller than DP on autosomes due to having large number of male samples, it seems reasonable to build VQSR models separately for X chromosome and autosomes. However, no such advice is proposed on GATK website.

**My question is**:

Should we build VQSR models separately for autosomes and X chromosome? If yes, we will have 4 VQSR models: Auto SNPs, Auto INDELs, X SNPs and X INDELs.

Cross-posting: https://www.biostars.org/p/9483110/

# Answer
You are correct that differences in ploidy between autosomes (always ploidy=2) and X chromosome mean that the same VQSR model cannot be applied to both cases equally.

However, this is the case only for male X chromosomes (in human species).

Therefore, yes, I'd recommend building a separate VQSR model when it comes to male X chromosomes. Importantly, the processing has to be different not only at the VQSR step, but from the very beginning at the HaplotypeCaller step (if you are using GATK also for variant calling), which has to be used with the option ploidy = 1.

So,

On one hand -\> Processing of autosomes (females and males) and x chromosomes (females only) can be done altogether with no differences all the way to VQSR (SNPs + INDELs as in best practices)

On the other hand -\> Processing of X chromosomes (only males) goes in separate, from variant calling to VQSR (SNPs + INDELs as in best practices)
> 2 votes

---
Tags: wgs, gatk
---


# Question
Title: Compressing the human genome to few megabytes
Multiple sources (see for instance this or this) discuss how genetic data will have scalability problems, given the huge file size of the human genome. The most straightforward encoding (see here) of the human genome requires about 700 Mb.

I came across this paper claiming to be able to store the human genome in about 4 Mb, having a reference genome and exploiting the fact that all human genomes are mostly equal. This paper is much older than the the other references discussing the scalability problems. Why is this technique not widely used?

# Answer
> 3 votes
*Answer from @alex-reynolds, converted from comments:*

As with any compression scheme, if you have to store a reference genome, along with the 4 Mb "diff", is it really compressed?

http://mattmahoney.net/dc/dce.html

For sequencing data, you might also look at the CRAM format, which also uses diffs against a reference genome to store reduced datasets. Labs are switching from BAM to CRAM to manage storage of increased numbers of sequencing datasets, with associated changes to pipelines and toolkits to use this format. For performance reasons, BAM may be preferred for analysis, and CRAM for long-term storage.

https://gatk.broadinstitute.org/hc/en-us/articles/360035890791-SAM-or-BAM-or-CRAM-Mapped-sequence-data-formats

https://www.ga4gh.org/news/cram-compression-for-genomics/

# Answer
> 1 votes
I saw this paper a few weeks back that claims to offer 4x more compression than typical formats for both fastqs and BAM/CRAM files.

Genozip - A Universal Extensible Genomic Data Compressor

https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/btab102/6135077

---
Tags: phylogenetics, human-genome
---


title: 10X scRNAseq: Sample mix-up

text:
# Question
Title: 10X scRNAseq: Sample mix-up
The student who was working on scRNA seq of KO and WT lines has made a mistake and he mixed both lines and generate the final sequencing data. Now, we are having gene expression data but don't know which is KO and which is WT.

It's a genetic KO. Looking at the presence of absence of the KO gene might not be helpful here as its not deep sequenced.

Cells were extracted from two stable line of animal; one wild type and other Chd2 gene knocked out. As a sequencing output I am having one fastq file in the format 28+10+10+90 . I can generate a count matrix out of it; but I dont know which barcode from WT and with one from KO.

Is there any way to segregate the output based on the KO gene expression?

Thanks

# Answer
If the samples were from different individual organisms you could try looking for mitochondrial variants that differentiate between them. You should have reasonable coverage of the entire mitochondrial mRNA transcriptome even in 10x datasets.

This tool might help: https://bioconductor.org/packages/devel/bioc/html/mitoClone2.html
> 1 votes

# Answer
*Answer from @haci, converted from comment:*

If your mutation of interest is "covered" by the "sample read", you should be able to assign which cell comes from which line. By "sampe read", I mean some 90+ bp read that is used for sequencing the insert (no the read for the cell barcode & UMI).
> 0 votes

---
Tags: scrnaseq, single-cell, gene-expression, 10x-genomics, cellranger
---


# Question
Title: Bash scripting FastQC for multiple fastq files in multiple directories
I am completely new to bioinformatics so I'm looking to learn how to do this.

I have multiple directories with fastq files: E.g; 10 Directories with each time series, each with Treatment and control directories, each with rep1 rep2 rep3.

For example: T9/Infected/Rep1/\*.fastq.gz.

I'm looking to create a loop to run fastQC on each fastq file instead of having to submit a separate job for each directory.

Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.

# Answer
`multiqc` kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the `multiqc` report and the actual fastqc `html` report to get a better idea of any error modes.

Going off of @Kubator's answer, I noticed that there was no command to run fastqc.

Here's a simple one-liner to run fastqc in parallel on all of your fastq files. The `-j 25` uses 25 threads. Change `25` to however many threads you want/have for max speed.

```
#Run fastqc on everything in parallel.
> find ../reads/ -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' | parallel -j 25 --verbose
# copies all the fastqc files to directory ./
> find ../reads/ -name '*fastqc.*' | xargs -I '{}' mv '{}' ./

```

These files might be output by `multiqc`, anyway!

edit - now the command is safe for whitespace in file names. Please, no newlines in filenames though!
> 5 votes

# Answer
This is the typical kind of job for snakemake.

Assuming your have one file per replicate named for instance `T9/Infected/Rep1/Rep1.fastq.gz`, you can prepare a file that you call `Snakefile` with the following content:

```
timepoints = list(range(10))
conditions = ["control", "infected"]
replicates = [1, 2, 3]

rule all:
    input:
        expand(
            "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html",
            time=timepoints,
            cond=conditions,
            rep=replicates)

rule do_fastqc:
    input:
        fastq = "T{time}/{cond}/Rep{rep}/Rep{rep}.fastq.gz"
    output:
        html = "T{time}/{cond}/Rep{rep}/Rep{rep}_fastqc.html"
    shell:
        """
        fastqc {input.fastq}
        """

```

Put this file in the directory that contains the `T*` directories and run `snakemake` from there.

The top `all` rule explains which files you want. The `do_fastqc` rule explains how to make one fastqc report from one fastq.gz file.

With a bit more work, this can be used to submit jobs to a computing cluster. Snakemake has some tools for this.

If you don't know the exact names of the fastq files but they all follow the same pattern, you will need to use the `glob` python module and do a little bit of programming to determine the possible values for `rep`, `cond` and `time`. The "snakefile" can contain any python code you want.

If there are no regular pattern in the file names, fix this issue first ;)
> 2 votes

# Answer
Example dir structure:

```
$ find FastQC/
FastQC/
FastQC/T9
FastQC/T9/Infected
FastQC/T9/Infected/Rep1
FastQC/T9/Infected/Rep1/test11.fastq.gz
FastQC/T9/Infected/Rep1/test1.fastq.gz
FastQC/T9/Infected/Rep2
FastQC/T9/Infected/Rep2/test2.fastq.gz
FastQC/T9/Infected/Rep3
FastQC/T9/Infected/Rep3/test3.fastq.gz

```

If understood well You need to run some job on every \*.fastq.gz file. Then You can do something like this (my example job is **gzip test**, replace with your job):

Rookie:

```
$ find FastQC/ -type f -name "*.fastq.gz" | xargs gzip -tv
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

```

Solid:

```
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} gzip -tv {}
FastQC/T9/Infected/Rep1/test11.fastq.gz:         OK
FastQC/T9/Infected/Rep1/test1.fastq.gz:  OK
FastQC/T9/Infected/Rep2/test2.fastq.gz:  OK
FastQC/T9/Infected/Rep3/test3.fastq.gz:  OK

```
1. find finds files with name \*.fastq.gz and outputs it with zero byte delimited (to support weird characters like space etc. in filenames)
2. xargs represents output as {} and passes it to gzip -tv

---

If You want to copy files inside one heap folder:

```
$ find FastQC/ -type f -name "*.fastq.gz" -print0 | xargs -0 -I {} cp -pv {} FastQC_heap/
`FastQC/T9/Infected/Rep1/test11.fastq.gz' -> `FastQC_heap/test11.fastq.gz'
`FastQC/T9/Infected/Rep1/test1.fastq.gz' -> `FastQC_heap/test1.fastq.gz'
`FastQC/T9/Infected/Rep2/test2.fastq.gz' -> `FastQC_heap/test2.fastq.gz'
`FastQC/T9/Infected/Rep3/test3.fastq.gz' -> `FastQC_heap/test3.fastq.gz'

```
> 1 votes

# Answer
You can use following command from the parent folder to run fastqc on `fastq.gz` files in subdirectories

```
find -name '*.gz' | xargs fastqc

```
> -1 votes

---
Tags: fastq, quality-control, bash, fastqc
---


# Question
Title: How to convert multiple single-end bam files to fastq using samtools
Hi I am trying to convert `bam` files generated from Ion Torrent Proton sequencing to `fastq` format so that I can upload them to KBase for analysis.

The files are named as follows: `01_thu159_IX_RNA_001.bam 01_thu163_IX_RNA_001.bam 01_thu164_IX_RNA_001.bam 01_thu190_IX_RNA_001.bam`

So I tried:

`samtools fastq *.bam > *.fastq`

However, that did not work. But, I know it is possible to use a placeholder for the name in code by using something like `$f` or I have also seen `%N%`. I am just not sure how to basically use those placeholders to be able to quickly apply a `for` loop, which will convert the files to fastq format without having to do each one individually. I have 28 samples each containing 4 `bam` files. Any advice on how I can do this would be greatly appreciated.

# Answer
> 2 votes
Using GNU parallel:

```
ls *.bam | parallel -j 2 "samtools fastq {} > {.}.fastq"

```

Use `-j` to control the number of parallel jobs. That only works for single-end data though as paired-end data must first be sorted or collated by name.

---
Tags: sam, samtools, fastq, format-conversion
---


# Question
Title: Snakemake rule issue with bcftools view filtering
I am creating a bioinformatics pipeline using Snakemake doing everying from adapter trimming of raw sequencing data in .fq to mapping, SNP calling and SNP filtering.

The issue i am having is when i try to filter specific chromosomes of the bcf files, using the --regions. The rule i have created below works and does create a filtered bcf file.

```
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} -o {output}"

```

However i would like to add the --regions parameter to the bcftools command. I have tried the shell part of the rule to several things

```
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions <chr1,chr2> -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {chr1,chr2} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1,chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions {{chr1} {chr2}}
 {{chr2}} -o {output}
bcftools view -Ob --threads {threads} -m2 -M2 -v snps {input.bcf} --regions ['chr1','chr2'] -o {output}

```

But it doesn't work, and i am having trouble finding any suggestions to this specific issue.

I know it works if i just run bcftools filter in my command line, so using --regions should work.

Could anyone have any suggestions how to add --regions to the snakemake rule?

# Answer
> 1 votes
Something like the following would work:

```
rule bcftools_filter:
    input:
        bcf = "/dir/{sample}_aligned_call.bcf",
    output:
        "/dir/{sample}_aligned_call_filter.bcf"
    params:
        regions="chr1,chr2"
    threads: 20
    shell:
        "bcftools view -Ob --threads {threads} -m2 -M2 --regions {params.regions} -v snps {input.bcf} -o {output}"

```

You could also include the regions as part of the output file name, rather than specifying it as a param.

---
Tags: snp, snakemake, bcftools
---


# Question
Title: Where to find rat genome for STAR alignment?
I'm trying to find the rat genome to carry out STAR alignment for several FASTQ files. Where would I find these files?

# Answer
> 2 votes
For example from Ensembl: https://www.ensembl.org/Rattus\_norvegicus/Info/Index

Can find other links via a search engine of choice.

---
Tags: sequence-alignment, genome
---


title: Assumptions of batch effect removal

text:
# Question
Title: Assumptions of batch effect removal
What does removing a batch effect (e.g. with `limma::removeBatchEffect`) assume about the batch effect?

Does it assume simply a constant batch effect for each level? E.g. if the levels are

```
Date1
Date2
Date3

```

Then one constant will be subtracted from all samples on Date3, another constant from all samples of Date2 and so on?

Or is it more complicated?

# Answer
You're basically subtracting a constant per-gene per-level. The relevant portion of the code is:

```
fit <- lmFit(x,cbind(design,X.batch),...)
beta <- fit$coefficients[,-(1:ncol(design)),drop=FALSE]
beta[is.na(beta)] <- 0
as.matrix(x) - beta %*% t(X.batch)

```

where `x` is your input matrix, `design` is the design matrix and `X.batch` is the matrix of batch and covariates. So you're fitting your data with the combined batch/design matrix (line 1), extracting out the batch coefficients (line 2), and subtracting out the expected effect (line 4). This is simple enough and remarkably effective for its intended use. You can get a bit fancier by using PCA or something else to instead allow individual samples to be more/less affected by a batch-effect, but added complexity has its own issues.
> 1 votes

---
Tags: limma, batch-effects
---


# Question
Title: Database for proteome-wide predictions of protein structures
Accuracies of protein structure predictors have improved quite a lot in recent years. Algorithms such as Rosetta have gotten robust enough to predict structures of large number of proteins. However, I could't find any initiative to make a database for proteome-wide predictions of protein structures or may be integrate the predicted structures in existing databases (PDB or Uniprot).

If any such exists, please let me know.

# Answer
> 4 votes
There are three initiatives I know of to have a go at this:

1. PConsFam, which collects data from this paper.
2. The Baker group's metagenomic study which you mention in the question.
3. The recent DMPfold work from our lab. This compares its results to the above two studies and discusses the effect on various model organisms.

The second and third of these don't provide a database per se, but the models are available for download and are linked to the Pfam database.

EDIT: things have moved on significantly since I first wrote this answer. As other answers mention, the AlphaFold Protein Structure Database is now the best resource in this area.

# Answer
> 4 votes
By far the most accurate model right now should be AlphaFold, and there are open structures across several proteomes available at https://alphafold.ebi.ac.uk/.

# Answer
> 3 votes
## Updated

As of 2021, the EBI alphaFold2 is the best resource. There are some caveats, which I have outlined in this blog post.

But majorly, it contains only monomeric predictions and not oligomers and complexes. As a result, for protein with a close homologue with a solved structure, it may be more handy to use a threaded structure as found in the SwissModel database (see original answer below) because it may be oligomeric and/or one can copy over ligands or binding partners if deemed suitable.

One can make one's own predictions with AlphaFold2, who have released a handy notebook. Also one can make complexes and oligomers with some tweaks or with the heavily improved notebooks found in ColabFold GitHub repo.

## Original Answer

It is worth saying that for model organisms, Expasy Swissmodel is the "classical" database of predicted proteome structures. It has several model organisms and gets updates regularly. The models are threaded only and the algorithm fell behind some time ago —their website however is modern and sleek. If using Rosetta afterwards, one should do a few cycles of FastRelax with a movemap preventing backbone movements or CA restrained, followed by a couple of cartesian bond/angle fixing cycles: the structure tends to blow up if energy minimised otherwise.

However, if one had a specific protein in mind predicting with RosettaCM, the I-Tasser server or Phyre server is by the recommended option. And actually, the Phyre server has a dataset of the precomputed human proteome.

About the PDB, the database has been purged of all *in* *silico* models and contains only empirical structures (NMR, X-ray, neutron diffraction and cryo-EM).

# Answer
> 1 votes
I am using DMPfold to predict protein structure of a 382 residue protein sequence. All goes wel, and even the files are produced but during the "seq2maps.csh", I get this error in between:

```
Running PSIPRED & SOLVPRED
Running PSI-BLAST with sequence c600m2-tail-before-lysin-sequence.temp.fasta ...
Segmentation fault (core dumped)
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.17 [Aug-26-2007]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for composition-based statistics:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingFATAL: Error whilst running blastpgp - script terminated!
Running PSICOV
Running FreeContact
Running CCMpred
 _____ _____ _____               _
|     |     |     |___ ___ ___ _| |
|   --|   --| | | | . |  _| -_| . |
|_____|_____|_|_|_|  _|_| |___|___|
                  |_|

using CPU (10 thread(s))

Reweighted 30404 sequences with threshold 0.8 to Beff=22917.9 weight mean=0.753779, min=0.0232558, max=1

Will optimize 62686429 32-bit variables

iter    eval    f(x)        ║x║         ║g║         step
1       1       1.50338e+07 152357      3.021629e+11    1.55e-06
2       1       1.47369e+07 152356      2.1447226e+11   1.18e-06
3       1       1.44292e+07 152355      1.4646483e+11   1.22e-06
4       1       1.41268e+07 152353      1.0388224e+11   1.5e-06
5       1       1.38295e+07 152358      7.337667e+10    1.92e-06

```

Can you please let me know if I need to rerun seq2maps.csh because of this error? And, how this error can be solved? Thank you!

---
Tags: database, proteins, protein-structure
---


# Question
Title: Error when using awk command to trim sequence files
I am attempting to edit my `fastq` sequence files with an `awk` command from a published article. The sequence files are all contained in one directory where I am running the `awk` command from the terminal in a Linux environment. There are 4 fastq files in the directory

```
01_thu159_IX_RNA_001.fastq
01_thu163_IX_RNA_001.fastq
01_thu164_IX_RNA_001.fastq
01_thu190_IX_RNA_001.fastq

```

I would like to trim all these files at once with the `awk` command so I wrote this `for` loop:

```
for filename in ./*.fastq*
do
awk -v s=10 -v e= 40 ' {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' ./${filename} >  ./${filename}_trimmed
done

```

However, I am getting the following error:

```
awk: fatal: cannot open file ` {if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;} ' for reading (No such file or directory)

```

I then thought maybe it was an error with the `for` loop so I tried the command as per the instructions in the article

```
awk -v s=10 -v e= 40 '{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' ./01_thu159_IX_RNA_001.fastq > ./Trim_01_thu159_IX_RNA_001.fastq

```

But this also produced the same error:

```
awk: fatal: cannot open file `{if (NR%2 == 0) print substr ($0, s+1, length ($0)-s-e); else print $0;}' for reading (No such file or directory)

```

I have checked the syntax for the `awk` command and it should be correct. Could anyone please give me advice as to what I am doing wrong and how to fix this issue?

# Answer
> 3 votes
You have a whitespace here:

```
e= 40

```

awk thinks that 40 is your rule and the rest of the line is a filename.

---
Tags: rna-seq, fastq, linux, awk, read-correction
---


# Question
Title: How do I include silent mutations in an oncoplot using maftools?
I have annotated files from annvar and I changed them into MAF by using `annovarToMaf` from ‘maftools’. Then I used `read.maf` to read them into R and then I used `merge_mafs` to integrate them together inside the last file.

I had several data frames but it separated the synonymous mutations (`maf.silent`) from the non-synonymous mutations, and I want to create an oncoplot using the `oncoplot` function for both synonymous mutations from the non-synonymous mutations.

But once I use this command I don't know how to integrate both data frames inside of the MAF file.

How can I create the oncoplot using `oncoplot` for both synonymous mutations and non-synonymous mutations?

# Answer
> 2 votes
This is late but I am putting it here for future reference. You can use the `vc_nonsyn` argument in the `read.maf` function to manually list the variant classifications you want to plot later. It will look like this:

```
df.maf <- read.maf(maf = df, vc_nonSyn = c("Frame_Shift_Del","Missense_Mutation","Nonstop_Mutation","Silent","Splice_Region"))

```

---
Tags: r, bioconductor
---


title: Alignment statistics

text:
# Question
Title: Alignment statistics
* Given a multi-fasta multiple sequence alignment file, is there a quick way/tool to calculate average pairwise sequence identity?
* Even better, given a whole genome alignment file (`.maf` or even `.hal`) is there a quick way/tool to calculate average pairwise sequence identity across the individual local multiple sequence alignments?

# Answer
There are tons of methods that do something similar to this. Is it necessary that it be *alignment* rather than a heuristic such as k-mer hashing?

If k-mer hashing is ok, then I would suggest mash or sourmash. There are additionally newer fancy techniques.

For just an old-school program that computes pairwise identities, PHYLIP has programs DNADIST and PROTDIST that exactly compute identities. Similar questions on biostars suggest that PHYLIP's utilities are recommended even relatively recently elsewhere. But it looks like Clustal Omega also computes identities for free.

For microbial genomes a common statistic is average nucleotide identity (ANI), which may be what you are looking for. There are a variety of tools to do this, here is one, and here is a webserver that will do it for you.
> 2 votes

---
Tags: sequence-alignment
---


# Question
Title: Pymol: color spheres with different colors
Hello I am a Pymol user

I want to color all alpha carbons (ca) and beta carbons (cb) from a certain ressidue in different color. How can I do this?

I selected the spheres with:

```
show spheres, C/1-9/ca
show spheres, C/1-9/cb

```

and if I use

```
set sphere_color, theColor

```

I color all the spheres what I don't want.

according to: `color gold, resi 200-‐221`
I tried `gold, spheres, C/1-9/ca` but this didn't work

# Answer
## Spheres

So how PyMOL stores how it shows atoms is that each atom has a colour and a represention, which is a binary number telling what to show. Therefore the colour for `spheres` is the same as the colour for `sticks`. The only way round this is to copy the object (`create copycat, original` and you set copycat differently)

To set the atom you do

```
select iliketheschaps, resi 200‐221 and name CA+CB
color tv_green, iliketheschaps
show spheres, iliketheschaps

```

The object `iliketheschaps` is a selection. Do make sure the range is a hyphen-minus (character next to 0) as opposed to an en-dash (–) which Word subs and on a mac you get with `cmd`+`-`

## non-PyMOL colours

Now you asked for "gold". This is not a colour in the menu —famous for the colour `limon`. Gold is a CSS colour name, but it won't work. You will have to colour by hex code. Say you want "gold" you need to google the colour name and hex. #FFD700.

```
color 0xFFD700, iliketheschaps

```

## The Fashion Police have been alerted

> This paragraph is in jest.

However, as you may know bioinformaticians are very fashion savvy for colours. Mustard was the in-colour from last year, so you really ought to use mustard (`0xffdb58`). This year is all about nectarine (`0xffa168`) (which is not salmon —this ain't 2015!). Coral (`0xff7f50`) and teal (`0x008080`) were from previous years and are still are going strong!
> 1 votes

# Answer
show sphere, A/1/C

set sphere\_color, red, A/1/C
> -1 votes

---
Tags: pymol
---


# Question
Title: polidy VCF file for Canvas somatic-wgs cnv calling
I am trying to run Illumina CANVAS cnv caller for Somatic-WGS. There is an option "--ploidy-vcf" which is mandatory to supply, but don't know what exactly that mean. I had supplied the CNV.vcf file from the germline WGS (matched normal run on Canvas smallPedigree), but it has given an error that says "Reference copy number \> 2 is not supported" at the last stage of run.

I would like to know your experience when you used CANVAS on somatic WGS samples. What ploidy-file is appropriate here? Should it be only chrX and Y in it? Your suggestions please.

# Answer
> 1 votes
It seems you need to supply vcf files specifying the ploidy of the sex chromosomes.

See https://github.com/Illumina/canvas/issues/89 for explanation and example vcf files.

---
Tags: wgs, illumina, cnv, somatic
---


title: Encoding different DNAs in pySBOL3?

text:
# Question
Title: Encoding different DNAs in pySBOL3?
I want to encode two different DNAs using pySBOL3. First dna is an engineered region usually composed by a transcriptional unit (TU) or a set of TUs. The second is vector which is a basically a plasmid that contains a set of the previously described dna. This is for differentiate between engineered on different plasmids.

Now for dna I used as type: `sbol3.SBO_DNA` and role: `sbol3.SO_ENGINEERED_REGION` and for vector I used as type: `sbol3.SBO_DNA` and role: `SO:0000155` (plasmid).

I would like to have feedback on this choice. For example the plasmid could also be an engineered plasmid `SO:0000637` or circular plasmid `SO:0002211` or just circular `SO:0000988`.

Here is my actual implementation:

```
for dna_id in params.get('dna'):
    #this needs to be handlded by vector
    dnas = Dna.object().filter(dna_id)
    dna = sbol3.Component(dnas.name(), sbol3.SBO_DNA)
    dna.roles = sbol3.SO_ENGINEERED_REGION
    dna.name = dnas.name
    dna.description = dnas.description
    if dnas == None:
        print('No dnas specified')

for vector_id in params.get('vector'):
    vectors =Vectors.object().filter(vector_id)
    vector = sbol3.Component(vectors.name(), sbol3.SBO_DNA)
    vector.roles = SO:0000155
    vector.name =vectors.name
    vector.description = vectors.description()

```

# Answer
Your choices match the same choices that I've made when designated types for inserts and vectors.

The more specific plasmid types are reasonable to use as well, but choosing just "plasmid" is certainly safe. That's why I have preferred to use that in my code because (1) my code might sometimes be applied to a different sort of plasmid, (2) I don't currently have a purpose for making the finer distinction, and (3) circular DNA would be indicated by SO.circular (SO:0000988) in any case.

*Side note: I would suggest using the tyto package so that your code can use names like `SO.plasmid` in your code instead of opaque symbols like `SO:0000155`.*
> 1 votes

---
Tags: phylogenetics, sbol
---


# Question
Title: Where to find publicly available Sanger chromatography data?
I am looking for any public available databases where I can download Sanger chromatography data ideally in Ab1 or SCF file format. I need extremely large amounts of this, since I want to use it for machine learning.

# Answer
> 2 votes
## NCBI Trace Archive

The NCBI Trace Archive is a permanent repository of DNA sequence chromatograms (traces), base calls, and quality estimates for single-pass reads from various large-scale sequencing projects. Contaning more than 2 billion traces.

---
Tags: ngs, sequence-analysis, public-databases, data-download
---


# Question
Title: Counting the number of gene isoforms in a GFF3, is this method correct?
Recently I've been tasked to count the numbers of gene isoforms for each locus in a `.gff3` file. I'm still doing my first steps in biology and bioinformatics, so I still struggle with some concepts and methodologies.

According to my understanding a gene isoform is a mRNA located in the same locus but with different start positions, so I did the following procedures:

1. Selected the entries in the `.gff3` that are classified as mRNA.
2. In my `.gff3` there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one.
3. Assigned each value from step 2. to their respective locus, resulting in the desired number of isoforms for each locus.

My question is *is this method correct*?

Why do I believe it could be wrong or inaccurate: I don't have a complete solid understanding of a gene isoform, and researching some more I've seen that not only the start point of the mRNA is relevant but other things such as the protein coding DNA sequence (which I believe could be identified in the file under the CDS tag) and untranslated regions (which I don't know if I can identify in my file), and in the method I described above I've only taken in consideration the Start position, what appears to be incomplete.

Additionally if the method is incorrect how can I take in consideration the additional circumstances, and if it's correct why didn't those circumstances influence the count so much?

# Answer
> In my .gff3 there is the additional information of the locus of the mRNA, so I checked for each locus how many unique values there were for the Start values and subtracted one;

Different isoforms of mRNA can share the same start coordinate. For example, they could share the same start coordinate but one skips an exon. MSpaint illustration below to help. Generally for GFF3 I would simply count e.g. how many mRNA have Parent=GENE0001234 to get the number of mRNA for GENE0001234, and you can expand that for all instances of GENENNNNNNNN

See also for example Wikipedia for Alternative splicing https://en.wikipedia.org/wiki/Alternative\_splicing
> 0 votes

# Answer
Gene isoforms are transcripts derived from the same gene that vary in their exons (blocks of the gene that make up the final transcript). This can alter the start position, end position, and/or internal sequence of the transcript. GTF/GFF files account for this by having multiple transcripts per gene (if applicable) that capture the various isoforms of the gene.

Here's a little example in R to count the number of transcripts for a few human genes found in an ENSEMBL reference GFF3 file.

```
library("GenomicFeatures")

gff3 <- makeTxDbFromGFF("http://ftp.ensembl.org/pub/release-104/gff3/homo_sapiens/Homo_sapiens.GRCh38.104.chr.gff3.gz")

> stack(vapply(transcriptsBy(gff3, "gene")[1:10], length, 1L))
   values             ind
1       5 ENSG00000000003
2       2 ENSG00000000005
3      16 ENSG00000000419
4       5 ENSG00000000457
5       9 ENSG00000000460
6       7 ENSG00000000938
7       6 ENSG00000000971
8       3 ENSG00000001036
9      14 ENSG00000001084
10      2 ENSG00000001167

```
> 0 votes

---
Tags: rna-seq, genomics, gff3
---