# Question
Title: Why I am getting different gene lengths for the same genes in different samples with rsem-calculate-expression?
I am running the `STAR-RSEM` pipeline, I have a simple shell script(\*), which I use to loop over fasta files corresponding to different samples. From the output, I am interested in the `*.genes.results` files, for which the columns are:

```
> names(my_df1)
[1] "gene_id" "transcript_id(s)" "length" "effective_length" "expected_count" "TPM" "FPKM"

```

To my surprise, whereas the `gene_id` and `transcript_id(s)` columns are the same, I see different `length` and `effective_length` values for the same genes in different sample results (omitting the `transcript_id(s)` column for formatting purposes):

```
> my_df1 <- fread(paste0("alignment/", my_files[1]))
> head(my_df1[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count   TPM
1: ENSG00000000003.14 2292.23          2242.23            978 26.69
2:  ENSG00000000005.5 1339.00          1289.00             11  0.52
3: ENSG00000000419.12 1024.45           974.45            469 29.45
4: ENSG00000000457.13 4204.00          4154.00            241  3.55
5: ENSG00000000460.16 2953.40          2903.40             62  1.31
>
> my_df2 <- fread(paste0("alignment/", my_files[2]))
> head(my_df2[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count    TPM
1: ENSG00000000003.14 2260.38          2210.38       10052.00 117.78
2:  ENSG00000000005.5 1269.22          1219.22         132.00   2.80
3: ENSG00000000419.12 1058.13          1008.13        4386.00 112.68
4: ENSG00000000457.13 4175.05          4125.05         873.00   5.48
5: ENSG00000000460.16 2783.92          2733.92         780.01   7.39
>
> my_df3 <- fread(paste0("alignment/", my_files[3]))
> head(my_df3[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count  TPM
1: ENSG00000000003.14 2662.59          2612.59             49 2.90
2:  ENSG00000000005.5  940.50           890.50              0 0.00
3: ENSG00000000419.12  956.87           906.87             58 9.89
4: ENSG00000000457.13 6364.00          6314.00             66 1.62
5: ENSG00000000460.16 2661.00          2611.00             15 0.89

```

The `length` and `effective_length` are defined per gene (and its transcripts) and should not differ across samples. Below is from the manual:

> A gene's 'length' and 'effective\_length' are defined as the weighted average of its transcripts' lengths and effective lengths (weighted by 'IsoPct')

I would appreciate any hints pointing out what I am missing.

* The script I am using:

```
for file in ...; \
do rsem-calculate-expression \
-p 8 \
--star \
# --output-genome-bam \
--star-path path-to-star-executable \
--star-gzipped-read-file \
--forward-prob=0 \
$file \
path-to-genome \ # should end with the string used for the "reference_name" parameter of "rsem-prepare-reference"
${file:11:7}; \ # sample prefix
done

```

EDIT: For those who will be interested, this entry from the RSEM Users Google Group demonstrates how to calculate gene length from weighted transcript length.

# Answer
If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.
> 2 votes

---
Tags: star, rsem
---


# Question
Title: Making some columns if a dataframe numeric
I have two dataframes

```
> head(a[1:4,1:4])
   Tumor_Sample_Barcode Chromosome Position End_Position
1:                   A1       chr4 90169866     90169866
2:                   A1      chr11 60235747     60235747
3:                   A1       chr1   983023       983023
4:                   A1       chr1 11346060     11346060
>

> head(c)
   Chromosome Position      VAF
1:       chrM     6691 0.610284
2:       chrM    14503 0.693325
3:       chr1 31412236 0.645161
4:       chr1 55693305 0.602941
5:       chr1 69963412 0.709302
6:       chr1 72608266 0.720000
>

> str(c)
Classes ‘data.table’ and 'data.frame':  44175 obs. of  3 variables:
 $ Chromosome: chr  "chrM" "chrM" "chr1" "chr1" ...
 $ Position  : num  6691 14503 31412236 55693305 69963412 ...
 $ VAF       : num  0.61 0.693 0.645 0.603 0.709 ...
 - attr(*, ".internal.selfref")=<externalptr>

```

I want to merge them like this but complaining about not being integer

```
> a=a[c, on = c("Chromosome","Position"), VAF := i.VAF]
Error in bmerge(i, x, leftcols, rightcols, roll, rollends, nomatch, mult,  :
  Incompatible join types: x.Position (character) and i.Position (integer)

```

I have tried something like these but not working

```
> c=c[!is.na(as.numeric(as.character(c$Position))), ]
>
> c=c[!is.na(as.numeric(as.character(c$VAF))), ]
>
> c$VAF <- as.numeric(as.character(c$VAF))

c$Position <- as.numeric(as.character(c$Position))

```

which not working

# Answer
> 2 votes
Your issue seems to be related to the column `Position` being a character column in `a`, you should try:

```
library(data.table)
a$Position <- as.numeric(a$Position)
a[c, on = c("Chromosome","Position"), VAF := i.VAF]

```

---
Tags: r
---


title: Seurat FindMarkers() output, percentage

text:
# Question
Title: Seurat FindMarkers() output, percentage
I have generated a list of canonical markers for cluster 0 using the following command:

```
cluster0_canonical <- FindMarkers(project, ident.1=0,  ident.2=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14), grouping.var = "status",
min.pct = 0.25, print.bar = FALSE)

```

And I got the following output, but I do not understand why the canonically expressed **ribosomal** markers seem to be expressed at high levels in both the cells within the cluster AND cells that are NOT in the cluster (i.e. both pct.1 and pct.2 are close to 1, but P value is highly significant). We would not expect this for canonical markers - we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster. I don't seem to have such an issue with non-ribosomal genes within cluster 0 and cluster 0 is the only cluster where I get a high number of ribosomal genes as canonical markers.

# Answer
Your interpretation of `pct.1` and `pct.2` is misguided. If `pct.1 == 1`, the given gene was detected in all of the cells. This carries no information about high expression or low expression in the cell, just that it was detected. `avg_logFC` is the indication of expression **levels**. So your expectation, that

> "we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster"

holds true, as indicated by the positive `avg_logFC` and significant p-values of the listed RP genes.

Seurat use nautral log, so the FC of RPS6 in cluster 0 vs. all other clusters indicated is 2.718281828459^.55947=1.750.

---

Seurat `FindMarkers()` documentation.
> 3 votes

---
Tags: seurat
---


title: Replace a prefix of a string

text:
# Question
Title: Replace a prefix of a string
I have a data frame like this

```
> head(a[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Position End_Position
1:                   A1       chr4 90169866     90169866
2:                   A1      chr11 60235747     60235747
    > dim(a)
[1] 753655    134
>
> unique(a$Tumor_Sample_Barcode)
 [1] "A1"  "A10" "A11" "A12" "A13" "A14" "A15" "A16" "A17" "A18" "A19" "A2"  "A20" "A21" "A22" "A23" "A24" "A25" "A26" "A27" "A28"
[22] "A29" "A3"  "A30" "A31" "A32" "A33" "A4"  "A5"  "A6"  "A7"  "A8"  "A9"

```

I want to change `a$Tumor_Sample_Barcode` from ***A1:A21*** to ***res1:res21*** but I don't how to do that

# Answer
I think you are looking for `gsub` function:

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

```

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

**Dummy data**

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

```

**Test of `gsub` vs `str_replace_all`**

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

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

# Answer
Other options:

* Drop 1st character and prefix with `"res"`
* We get 5x speed up if we set `fixed = TRUE` for *gsub*.

```
microbenchmark::microbenchmark(
  # dc37 answer
  gsub("A","res", dt$Barcode),
  stringr::str_replace_all(dt$Barcode, "A", "res"),
  # other options
  paste0("res", substring(dt$Barcode, 2)),
  gsub("A","res",dt$Barcode, fixed = TRUE),

  # bench set up
  check = "identical",
  unit = "relative"
  )

Unit: relative
                                             expr      min       lq     mean   median       uq      max neval
                     gsub("A", "res", dt$Barcode) 5.171682 5.132816 4.283014 5.033968 3.177772 2.545173   100
 stringr::str_replace_all(dt$Barcode, "A", "res") 6.914598 6.870849 5.747396 6.728126 4.809257 3.518024   100
          paste0("res", substring(dt$Barcode, 2)) 1.983757 1.991876 1.779557 1.968488 1.480892 3.415619   100
       gsub("A", "res", dt$Barcode, fixed = TRUE) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000   100

```
> 7 votes

---
Tags: r
---


title: Read-informed statistical phasing

text:
# Question
Title: Read-informed statistical phasing
I'm wondering if there are any tools out there that perform statistical phasing of genetic data (e.g. Eagle) that take read-based phasing from GATK into account. I've found this paper which sounds like what I am looking for, but that tool was developed specifically for Hi-C data, and I only have exomes.

Thanks

# Answer
shapeit4 is currently (probably) the most accurate and quickest phasing algorithm which can account for sequencing reads in the model. They suggest you extract the phase information from the reads using WhatsHap first however. The guide for WhatsHap can be found here and has some info on using data processed using GATK.
> 1 votes

---
Tags: exome, genetics, phasing
---


# Question
Title: Module preservation and hub genes finding
I ran module preservation analysis between two condition and on the basis of preservation z summary i would to go for downstream analysis.So I m interested in looking for hub-genes of those modules which have lowest z summary. But so far in the tutorial i didn't find a section in the WGCNA tutorial where how to find hubgenes after module preservation analysis was mentioned.

Any suggestion would be really appreciated

# Answer
> 1 votes
Simply select the module(s) you are interested in and look for hub genes using the standard calculations. Having done a module preservation analysis does not change the procedure in the slightest.

---
Tags: networks, wgcna
---


# Question
Title: Any user friendly way to find rare mutations in whole genome raw?
Is there any user friendly way to find rare mutations in the individual human whole genome sequencing raw data? (from Dante, 30x coverage).

To be more specific, I want to find mutations from this paper:

https://docs.google.com/document/d/1EkRMuD6J0-zyMY3MegKhz7Th0hzhzVXg6DWg2Hm6l2Y/edit?usp=sharing (it's very short, less than one page).

But I'm confused. In their paper they just talk about genes (RAD21, B3GAT2, SMC3, SCN11A , SCN5A, SCN9A, SCN10A, SCN11A, TRPA1), but not mutations? Can we find that diseased-mutations which they talking about in genes, which their list? Or is there not enough data in the paper/study for this?

And if there is not enough data, which data I need to request from authors?

Or those genes have kind a "gold standard sequences" and if that sequence different from standard - there is "diseased" gene? (I have only very basic genetics and bioinformatics knowledges)

# Answer
You need to run a SNP finding pipeline, filter for SNPs in those genes and then run a SNP annotation pipeline. (You could annotate before you filter but it will take longer).

SNP finding usually has three steps, 1) read alignment to reference, 2) SNP calling and 3) SNP filtering (by quality).

You then bring in gene annotation data to filter SNPs to just those genes.

You then run SNP annotation using tools like VEP, SIFT, PolyPhen, etc.

People use tools like GATK to do this.

Hopefully that's enough information to get you started.
> 1 votes

---
Tags: gene, ngs, genome, sequence-analysis
---


# Question
Title: Genetic Genealogy related to SNP variants and INDEL variants
**Background:** I am looking for some advice on the use of SNP variants versus INDEL variants for determining Genetic Genealogy, specifically for determining accurately the "START" and "END" of shared DNA segments (shared between two DNA testers' results). The test data used is RAW DNA from direct to consumer DNA testing companies such as AncestryDNA and 23andMe.

**Question:** For the purpose of Genetic Genealogy, are the INDEL variants as useful as the SNP variants? Should the INDEL variants be ignored? Or, alternatively, do the INDEL variants hold a greater significance than the SNP variants?

I hope that using INDEL variants would increase the data in the comparison algorithm, thus increasing the certainty and resolution when determining each shared segment.

# Answer
InDels are just as informative, but are less reliable to call. Errors should be averaged out in this kind of analysis though (so long as they don't predominate).

Both SNPs and InDels are valid data points to consider when comparing triads, but identifying IBD segments will only ever be probabilistic.
> 1 votes

# Answer
first of all there's no way to accurately determine the start/end position of shared DNA segments when using microarray data from DTC's like 23andMe and Ancestry.

INDEL's are indeed important when it comes to identify shared segments though I'm not sure how often the change over the course of several generations.

My mum's sister has an INDEL that both my mum and her brother didn't inherit, it will be inherited by her children (my 1st cousins).

Unfortunately the microarray from 23andMe doesn't give you any detailed information of than DI, II and II but not how long the deletion or insert is.

To determine such it requires to use WGS data.
> 0 votes

---
Tags: genetics
---


# Question
Title: Tool for calibrating base quality scores?
A given sequencing machine assigns a 'certainty' to each base call in the form of a per-base quality score (FASTQ format).

I have reads from several machines aligned to the corresponding reference (BAM format or similar). Ignoring SNPs and misalignments (for simplicity), mismatches should represent sequencing errors.

I'd like to parse through the FASTQ and BAM file and evaluate how closely the per-base quality score reflects 'true' sequencing errors.

Is there an existing efficient, user-friendly, open source, well documented tool to do this? Failing that, is there a paper doing this that claims to have made the tool freely available?

Another method would be to look at bases corrected using kmer analysis or similar. Again, the per-base quality score can be compared to the 'true' sequencing errors.

Any tool to do that? I guess it could be baked into an existing error correction tool.

Do short read aligners actually make use of per-base quality score?

Many thanks,

# Answer
> 1 votes
I think what you are looking for is GATK's Base Quality Score Recalibration. Their documentation is great so I won't rewrite it here.

> BQSR stands for Base Quality Score Recalibration. In a nutshell, it is a data pre-processing step that detects systematic errors made by the sequencing machine when it estimates the accuracy of each base call.

---
Tags: sam, fastq, quality-control
---


# Question
Title: Processed spliced and unspliced count matrices for existing scRNA-seq atlases
Before I do this myself...

Is there any compiled batch-corrected dataset of the major scRNA-seq atlases (Mouse Cell Atlas, Mouse Organogenesis Atlas, Mouse Gastrulation Atlas, Tabula Muris, Tabula Muris Senis, etc.) with processed spliced AND unspliced count matrices?

I would be surprised if there is, but in the event there isn't, are there new pipelines in development that will make processing of raw data easier for studying RNA velocity?

# Answer
> 1 votes
I am not sure if there are preprocessed datasets available, however, Kallisto + Bustools was designed specifically for the rapid quantification of spliced and unspliced transcripts.

---
Tags: r, python, scrnaseq
---


# Question
Title: 'Wildcards' object has no attribute 'sample'
I'm using Snakemake and I'm trying to get all my QC reports into multiqc but I get the following error:

```
WorkflowError in line 120 of /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full/Snakefile:
'Wildcards' object has no attribute 'sample'

```

Rule 120 is the start of my multiqc rule which looks like this:

```
rule multiqc:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.stat', '{sample}.cnt'),
                os.path.join(rseqc_dir, '{sample}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{sample}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{sample}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{sample}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{sample}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.junction.xls'),
                os.path.join(rseqc_dir, '{sample}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{sample}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{sample}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{sample}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{sample}.GC.xls'),
                os.path.join(rseqc_dir, '{sample}.NVC.xls'),
                os.path.join(rseqc_dir, '{sample}.qual.r'),
                os.path.join(rseqc_dir, '{sample}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{sample}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastqc_dir, '{sample}_R{read}_001_fastqc.html')],
                sample=samples['samples'], read=reads)
    output:
        os.path.join(qc_dir, 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        multiqc . --outdir qc --ignore .snakemake --force
        '''

```

My "rule all looks like this:"

```
rule all:
    input:
        expand([os.path.join(analyzed_dir, '{sample}.genes.results'),
                os.path.join(qc_dir, 'multiqc_report.html')],
                sample=samples['samples'])

```

I don't see anything wrong with it, in fact, it is very similar to my "rule all" I used before and which worked. Can anyone see an error?

EDIT

I limited the rule to only one input, I still have the same error:

```
rule multiqc:
    input:
        expand(os.path.join(rseqc_dir, '{sample}.geneBodyCoverage.txt'),
                sample=samples['samples'])
    output:
        os.path.join(qc_dir, 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        multiqc . --outdir qc --ignore .snakemake --force
        '''

```

This is the ouput, you can see it read the samples correctly:

```
run_snakemake.sh
/rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake/lib/python3.6/site-packages/snakemake/workflow.py:12: FutureWarning: read_table is deprecated, use read_csv instead, passing sep='\t'.
  from collections import OrderedDict
Building DAG of jobs...
/rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake/lib/python3.6/site-packages/snakemake/workflow.py:12: FutureWarning: read_table is deprecated, use read_csv instead, passing sep='\t'.
  from collections import OrderedDict
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 256
Job counts:
    count   jobs
    1   all
    1   multiqc
    2

[Mon Jan 14 17:46:40 2019]
rule multiqc:
    input: qc/rseqc/0067_P2018SEQE27R01_S1.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R02_S2.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R03_S3.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R04_S4.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R05_S5.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R06_S6.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt, qc/rseqc/0067_P2018SEQE27R08_S8.geneBodyCoverage.txt
    output: qc/multiqc_report.html
    jobid: 2

WorkflowError in line 119 of /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full/Snakefile:
'Wildcards' object has no attribute 'sample'
17:46 nlv24077@kiato /rst1/2017-0205_illuminaseq/scratch/swo-406/test_snakemake_full > ll qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt
-rw-rw-r-- 1 nlv24077 2017-0205_illuminaseq 303 Jan  7 16:00 qc/rseqc/0067_P2018SEQE27R07_S7.geneBodyCoverage.txt

```

EDIT 2:

This is how I start Snakemake:

```
> cat run_snakemake.sh
source activate /rst1/2017-0205_illuminaseq/scratch/swo-406/snakemake
snakemake --dag | dot -T svg > dag.svg
snakemake --cluster-config cluster.json --cluster "qsub -l nodes={cluster.nodes}:ppn={cluster.ppn} -N {wildcards.sample}_{rule}" --jobs 256 --use-conda --rerun-incomplete

```

EDIT 3:

Ok, that last remark (by @finswimmer) triggered me, I removed the "{wildcards.sample}\_" from the above command and then it works. Or, works, I get a very different error, it is multiqc specific and solved as specified here: https://multiqc.info/docs/#locale-error-messages. When adding the two lines as suggested this error disappears and everything works as expected, even with the long inputs list.

```
RuntimeError: Click will abort further execution because Python 3 was configured to use ASCII as encoding for the environment. Consult https://click.palletsprojects.com/en/7.x/python3/ for mitigation steps.

```

Anyway, I'm guessing Snakemake does not understand how to name this task when it uses an entire list of {sample}'s... Is this a bug? What should be used? I guess it is impossible to use {sample} for naming for this rule. It's a shame I have to kill the naming for all rules because using sample\_rule is very nice for keeping track.

# Answer
Snakemake gets `wildcards` from the parsing the input/output file names. For example, if you had

```
rule sort:
    input: "{file}.bed"
    output: "{file}.sorted.bed"
    shell: "sort -k1,1 -k2,2n {input} > {output}"

```

then the `wildcards` variable would be set as `wildcards = {"file": ...}`.

In your `multiqc` rule, you have the `{sample}` variable within your `expand()`, but this variable doesn't exist outside of `expand`'s local context. You don't actually have any wildcards there, nor in your output file (`qc_dir` and `rseqc_dir` are variables you define elsewhere, but not "wildcards" as Snakemake interprets them).

So you're getting the error `'Wildcards' object has no attribute 'sample'` because there are no wildcards for that rule.

Your rule is currently set up to process every sample that you have listed in `samples['sample']`. If you want your `multiqc` rule to work for a single sample, simply remove the `sample=samples['sample']` and double brace the `{sample}` wildcard.

```
rule multiqc:
    input:
        expand([os.path.join(analyzed_dir, '{{sample}}.stat', '{{sample}}.cnt'),
                os.path.join(rseqc_dir, '{{sample}}.bam_stat.txt'),
                os.path.join(rseqc_dir, '{{sample}}.clipping_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.deletion_profile.txt'),
                os.path.join(rseqc_dir, '{{sample}}.infer_experiment.txt'),
                os.path.join(rseqc_dir, '{{sample}}.geneBodyCoverage.txt'),
                os.path.join(rseqc_dir, '{{sample}}.inner_distance.txt'),
                os.path.join(rseqc_dir, '{{sample}}.insertion_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.junction.xls'),
                os.path.join(rseqc_dir, '{{sample}}.junctionSaturation_plot.r'),
                os.path.join(rseqc_dir, '{{sample}}.mismatch_profile.xls'),
                os.path.join(rseqc_dir, '{{sample}}.read_distribution.txt'),
                os.path.join(rseqc_dir, '{{sample}}.pos.DupRate.xls'),
                os.path.join(rseqc_dir, '{{sample}}.seq.DupRate.xls'),
                os.path.join(rseqc_dir, '{{sample}}.GC.xls'),
                os.path.join(rseqc_dir, '{{sample}}.NVC.xls'),
                os.path.join(rseqc_dir, '{{sample}}.qual.r'),
                os.path.join(rseqc_dir, '{{sample}}.RNA_fragment_size.txt'),
                os.path.join(rseqc_dir, '{{sample}}.STAR.genome.sorted.summary.txt'),
                os.path.join(fastqc_dir, '{{sample}}_R{read}_001_fastqc.html')],
                read=reads)
    output:
        # this {sample} is in single {} and is recognized as a wildcard
        os.path.join(qc_dir, '{sample}', 'multiqc_report.html')
    conda:
        "envs/multiqc.yaml"
    shell:
        '''
        #!/bin/bash
        # using {wildcards.sample} here since MultiQC likes directories
        multiqc . --outdir qc/{wildcards.sample} --ignore .snakemake --force
        '''

```

When wildcards are recognized, you'll see them in the job list.

Running `snakemake 0067_P2018SEQE27R01_S1/multiqc_report.html`, for example, should give you something like:

```
rule multiqc:
    input: qc/rseqc/0067_P2018SEQE27R01_S1.geneBodyCoverage.txt, ...
    output: qc/0067_P2018SEQE27R01_S1/multiqc_report.html
    jobid: 2
    wildcards: sample=0067_P2018SEQE27R01_S1

```
> 5 votes

# Answer
Although @JamesHawley answer was informative, I don't think it solved the initial issue requested by the OP. (**edit**: or maybe it partially did, I'm not sure as I got a bit confused with the wildcards on the initial question/answer. I'm anyway leaving this answer here, in case another way of explaining it is useful for others)

### Solution

I think you only need to change the output of your `rule multiqc` from:

`os.path.join(qc_dir, 'multiqc_report.html')`

to

`os.path.join(qc_dir, '{sample}.html')`

This seemed counter-intuitive to me initially, but I think what's happening is that `snakemake` will look at `rule all` and see that you want an output called "qc/multiqc\_report.html". And the only place where it can find something of that sort is in `rule multiqc`, by replacing the wildcard `{sample}` by the string `multiqc_report`.

It felt confusing to me, because the same `{sample}` wildcard is used in other rules (namely the rules you use to generate the outputs in your `analysed_dir`). But `snakemake` will be smart to replace `{sample}` with whatever it needs to recreate the output files you're requesting.

---

### Longer explanation

Here is a (hopefully) reproducible example, if it helps.

First create some fake directory structure with some files:

```
mkdir inputs
mkdir analysed_dir
mkdir qc_dir

touch inputs/file1.txt
touch inputs/file2.txt

```

Then, this would be the snakefile:

```
# Define your variables, fetch things from config, or csv, etc...
input_dir = "inputs"
output_dir = "outputs"
samples = ["file1.txt", "file2.txt"]

# Rule all defining all the output files required
rule all:
  input:
    expand("analysed_dir/{samples}", samples = samples),
    "qc_dir/multi_qc_report.html"

# this rule is equivalent to whichever rules process your data individually
rule analyse_each_file:
  input:
    "inputs/{sample}.txt"
  output:
    "analysed_dir/{sample}.txt"
  shell:
    "cat {input} > {output}"

# This rule is equivalent to your multiqc step where multiple files result in a single output
rule multiqc:
  input:
    expand("analysed_dir/{samples}", samples = samples)
  output:
    "qc_dir/{sample}.html"
  shell:
    "cat {input} > {output}"

```

to generate this pipeline:

And then issuing:

`snakemake --cluster "qsub -N {wildcards.sample}-{rule}" --jobs 10`

Generates the following job names:

* file1-analyse\_each\_file
* file2-analyse\_each\_file
* multi\_qc\_report-multiqc

**Edit with further explanation:**

Notice in the DAG how the `sample` wildcard takes different values for the different rules. I think the snakemake logic is as follows:

* Starting with **`rule all`** I see that you want three files as input (the first two come from the `expand()` call):
    + `analysed_dir/file1.txt`
    + `analysed_dir/file2.txt`
    + `qc_dir/multi_qc_report.html`
* Let me start with the first one `analysed_dir/file1.txt`. The file doesn't exist, so let me look in the other rules and see if I can find it in their output. Actually in **`rule analyse_each_file`** I see the output is set to `"analysed_dir/{sample}.txt"`, so if I replace `{sample}` by `file1.txt` I can recreate your file from it. Therefore, the input for that rule will be `"inputs/file1.txt"` (and this file already exists, so I can issue this job).
* For the second file `analysed_dir/file2.txt`, well it's a similar case as with the previous one, so I'll issue a job with the same **`rule analyse_each_file`**, but this time `{sample}` = `file2.txt`.
* Finally, you want a file called `qc_dir/multi_qc_report.html`. Where can I create such a file from? Well, the output of **`rule multiqc`** is set to `"qc_dir/{sample}.html"`, so if I make `{sample}` = `multi_qc_report`, then I can issue the job and produce your file of interest. In this case, the rule takes as input two files `"analysed_dir/file1.txt"` and `"analysed_dir/file2.txt"`, which I already know how to create from another rule.

The key thing to realise is that `snakemake` **defines the wildcard `{sample}` for each rule separately**. So, in this example, `wildcard.sample` in `rule analyse_each_file` does not take on the same values as in `rule multiqc`.
> 0 votes

---
Tags: snakemake
---


title: Genome annotation of a Bacillus strain

text:
# Question
Title: Genome annotation of a Bacillus strain
I finished the sequencing of the genome of a biocontrol activity strain *Bacillus amyloliquefaciens* by Illumina hiseq. I have already gone through some papers but there were no details of processing in it.

Is anywhere I can find some resources to get my research going on to genome annotation, compare it with other microbial genome data, and analyze non-ribosomal peptide synthetase, find out more genomic information about biocontrol activity?

For example, how to locate 16S DNA and some lipopeptide antibiotic genes (fungycin gene, surfactin gene) in my bacterial genome, and compare it with the reference strain and draw a phylogenetic tree. And how to use my bacterial genome and 2 reference genomes to draw a genomic circle diagram.

# Answer
> 2 votes
You said that you want to compare with other microbial genome data. Have you already annotated your genome? If not, you might consider a pipeline I've developed, I can gladly help if you face some problems using it.

It is called bacannot, it does a complete annotation using Prokka, Kofam, CARD, VFDB, and other resources. It might help you to have a first scratch of your data since it also collapses annotation information in a genome browser

---
Tags: sequence-annotation, genome, sequence-analysis
---


# Question
Title: How to generate two fasta files from samtools phase output?
samtools phase is an easy-to-use tool to phase variants from a bam file and a reference.

For a sequence that has two haplotypes, like this:

```
   position: 123456789
haplotype 0: GACCATATT
       SNPs:  **    *
haplotype 1: GTACATACT

```

The output looks like this:

```
CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
CC
CC
PS      frame-1   2    9
M1      frame-1  60    2      A       T       1       3       0       23      0
M1      frame-1  60    3      C       A       2       35      1       24      0
M1      frame-1  60    8      T       C       3       46      0       30      0

```

How can I take the reference fasta file and the output of `samtools phase` to generate one `fasta` file for each of the haplotypes?

# Answer
**Answer January 20th, 2020.** If you run `samtools phase` with the option `-b my_prefix` it outputs a bam file for both haplotypes, named `my_prefix.0.bam` and `my_prefix.1.bam`. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the **January 14th, 2020** answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that `samtools phase` found. `samtools phase` does not find indels.

**Answer January 14th, 2020.** This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.

Run the program like this:

`samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta`.

Here is the code of `phase2fasta.py`:

```
#!/usr/bin/env python

import os
import sys

if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
    raise Exception("The first and only argument is the reference fasta.")
allele0 = []
seqname = ""
with open(sys.argv[1]) as f:
    for line in f:
        if line[0] == ">":
            seqname = line[1:].split(" ")[0].strip()
        else:
            for c in line.strip().upper():
                allele0.append(c)
allele1 = [x for x in allele0]
for line in sys.stdin:
    if line[0:2] == "M1":
        fields = line.strip().split()
        if fields[1] == seqname:
            pos = int(fields[3]) - 1
            allele0[pos] = fields[4]
            allele1[pos] = fields[5]
print(">{}_allele0".format(seqname))
print("".join(allele0))
print(">{}_allele1".format(seqname))
print("".join(allele1))

```
> 3 votes

---
Tags: sam, fasta, samtools, phase
---


# Question
Title: Seurat define manually a cluster and find markers
I am aware of this question Manually define clusters in Seurat and determine marker genes that is similar but I couldn't make tit work for my use case.

So I have a single cell experiments and the clustering id not great I have a small groups of 6 cells (I know it is extremely small, but nonetheless I would like to make the most of it) that are clearly isolate in UMAP and display marker that I am interresting in.

I managed to get those cells names in barcode (e.g.: AGCCAGCTCGCTTTAT).

I would like to compute differential expression of those cells against all others groups (using FindMarkers or FindAllMarkers).

I tried to manually define the cluster in: Seurat\_obj@meta.data$seurat\_clusters But using FindMarkers (it says this groups don't exist)

then I try to modify s\_19270BC@active.ident

but then it could not found metadata attached to my newly define group.

I am kind of stuck and any help is warmly welcome

# Code

My R is quite rusty...

```
names_ # a vector with the names of the cell I want as a cluster

# my seurat object have 5 cluster 0..4. So I want to add a 6th clusters with id 5

#First attempt manually modify the vector active.ident and cluster
# they are both factor so we need to add a factor level
levels(seurat_obj@active.ident) = c(levels(seurat_obj@active.ident), 5)
levels(seurat_obj@meta.data$seurat_clusters) = c(levels(seurat_obj@meta.data$seurat_clusters), 5)
# then modify them :
seurat_obj@active.ident[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5
seurat_obj@meta.data$seurat_clusters[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5

FindMarkers(s_19270BC, ident.1 = 5)

```

So This morning it works... I don't know what I have done yesterday...

EDIT: I added the Code I used

# Answer
> 2 votes
It would have been better if you have shared your code as well as the error messages you got format these as such.

As far as I can understand, your problem is that you failed to make your cluster as default ids. For this you can use:

```
seurat_obj <- SetIdent(seurat_obj, id=seurat_obj@meta.data$column_of_interest).

```

The cluster identities are stored in the `@meta.data` slot of your object.

---
Tags: r, seurat, differential-expression
---


# Question
Title: Appending some columns from one tabular file to another
I have a list of files in which I have this information

```
Chrom   Position        Ref Alt TumorVAF    NormalVAF
chr1    56987           CA  C   0.222222    0
chr1    4133415         AG  A   0.425       0
chr1    4611030         A   AC  0.36        0

```

This is the name of one of my files

```
LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf.parsed.txt

```

I have another matched file which has too much information

This is the name of these big files

```
LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf_fixed_vcf.txt.hg19_multianno.txt

```

As you are seeing, most of the name of txt files are matched

I want to cut the 4th and 5th columns from my small files and add these columns to the matched big files

Manually this is possible but time consuming for a lot of files

Is these any way to do that by coding?

# Answer
> As you are seeing, most of the name of txt files are matched

Our eyes and brains are smarter than computers. We can see patterns that can be difficult to describe. For a computer to be useful in this task, the difference between the two file names must be explicitly defined.

In your example (which seems to be an ANNOVAR-generated output file with a prefix specified during annovar run), the common part seems to be the output of `$(basename LP6005334-DNA_H01_vs_LP6005333-DNA_H01.passed.somatic.indels.vcf.parsed.txt .parsed.txt)`. If this is true for all files, you could automate that using a loop such as:

```
for f in *.parsed.txt
do
    other_file=$(basename "$f" .parsed.txt)
    # Rest of your logic goes here
done

```

You could then `paste` the columns from `$f` to `$other_file` like so (not tested):

```
paste "$other_file" <(cut -f3,4 $f)

```

This is an outline to get to your solution. You should be able to build off this.
> 2 votes

---
Tags: r, bash, loop
---


# Question
Title: Ribosomal RNA QC quantification using SILVA
I would like to create a simple QC check of RNA-Seq data that simply maps the data to rRNA and then counts the number of reads that map.

I've done this manually for human using sequences for 5S, 5.8S, 18S, and 28S rRNA (rDNA, actually). I gathered those sequences manually to create a small fasta file to use as a reference and it worked well.

I would like to have a more definitive source for this reference though, and I don't necessarily want to use a species-specific database, so I thought that using a fasta file downloaded from SILVA would be great, but I don't know if I'm using the wrong fasta file or what. I'm not very familiar with SILVA, but I noticed that when I tested it using the same sequencing data, I got unexpectedly many fewer hits than I did with my manually created 4-sequence reference than with the SILVA NR99 fasta file.

So I blasted my 4 rDNA sequences against the SILVA NR99 data and only the 18S was found. The others (5S, 5.8S, and 28S) don't appear to be there. Am I misunderstanding the nature of the SILVA project or could I be doing something wrong?

*UPDATE*: Please note, this question is about **where to find rRNA and/or rDNA sequences to use as a reference**, not about tools that can be employed in quantifying rRNA read content. While I appreciate hearing about tools that can generate nice reports about rRNA reads in an RNA-Seq run, the question is not answered unless it addresses the missing rRNA/rDNA sequences that are necessary for comprehensive quantification. Any tool can count rRNA/rDNA mappings, but the results will only be as good as the rRNA reference data supplied (unless it is a reference-less tool). And it's the lack of 5S, 18S, and 28S ribosomal subunits in SILVA that is my main concern. It's the difference between 1k reads and 500k reads mapping to rRNA/rDNA sequences.

# Answer
> 1 votes
SILVA is intended for metagenomic samples and thus does not contain the human ribosomal subunits. While the 5.8S subunit is in the ensembl, ncbi, and ucsc annotations, and there are numerous rRNA "biotype" (or "biomol") entries, the 5S, 18S, and 28S subunit sequences are absent. Even downloading rRNA entries in the rmsk table from ucsc will not return all 4 major human subunits.

You can however (as I discovered) download all rDNA gene sequences (including those 4 subunits) from the nucleotide db on ncbi using the search term "biomol\_rrna\[PROP\]". This returns 240,240 coding sequences (or 30,775 gene sequences). You can then use this as a reference for mapping, blasting, or as input to rRNA QC tools.

The following e-Utilities commands from ncbi will download the gene sequences:

```
esearch -db nuccore -query "biomol_rRNA[prop]" | efetch -format fasta > all_rrnas.fasta

```

This is not a common analysis. Usually, rRNA reads can just be ignored and analyses are unaffected, which is likely why there is a paucity of records for these sequences in annotation files, but mapping to just these sequences could make for a simple QC measure as a gauge for rRNA contamination. While the sequences may be present in the genome, they're not in the annotations and may even be masked as repetitive, since there are so many copies. So you have to be careful to count unique reads and not mappings.

# Answer
> 0 votes
I'd recommend using kraken2 for this. There's a pre-built kraken2 database for SILVA (16S only) which can be downloaded from here:

https://ccb.jhu.edu/software/kraken2/downloads.shtml

For other ribosomal subunits, it's necessary to create a kraken2 database from scratch. 16S/18S subunits are in the SSU database, while 23S/28S are in the LSU database:

https://www.arb-silva.de/download/arb-files/

Note that SILVA is a database of collected rRNA sequences. As far as I'm aware there's no canonical "human" 18S or 28S sequence in SILVA, but it includes a collection of multiple 18S & 28S sequences that have been annotated with the human taxon.

Kraken2 has a command-line interface for creating reports:

```
kraken2 --db <location_of_database> --gzip-compressed input.fastq.gz --report output.kraken.report.txt > output.kraken.txt

```

Pavian can be used for visualising the results and creating proportion tables.

# Answer
> 0 votes
SortMeRNA has fasta files from RFAM and Silva for 5S, 5.8S, 18S, and 28S. You can download them from GitHub

---
Tags: ribosomal
---


# Question
Title: Mapping a list of cells in seurat featureplot
I have 209 cells, I clustered them by Seurat to 4 clusters. By Featureplot I am able to track a gene in clusters:

Higher color shows higher expression.

Now, for some genes I want to highlight some cells in Featureplot so that apart from yellow or red colours I want to colour a subsets of cells with another color. I mean I want to map a list of cells in Featureplot or tSNE plot. Let's say I want to know the location of cells 1, 4, 80 and highlight them with another color. My Seurat object in this link.

Seurat itself beautifully maps the cells in Featureplot for defined genes with a gradient of colours showing the level of expression. Saying I have genes A and B, in excel.

I have coloured cells that express a gene \> mean + se, \< mean - se or between these values. For instance, for this gene, 36 cells express this gene \> mean + se, I want to map these cells in Featureplot or tSNE plot in distinct colour so I can locate them in clusters easily. Something like binary (on off) expression to relative expression.

How I can change the colour of some cells based on my threshold for this gene please?

# Answer
> 5 votes
To color the TSNEPlot, you can generate a new column in metadata with the expression levels (High, low, etc). Then use pt.shape to set a shape for each identity.

To show binary expression based on expression you first have to define the list of cells that are below or over your threshold. Once you have those lists you can use SetIdent() in Seurat to color those groups.

```
load(file= "~/seuset_16.RData")

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
seuset_h16 <- StashIdent(seuset_h16, save.name = "expr")
seuset_h16 <- SetAllIdent(seuset_h16, id="res.1")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, pt.shape ="expr")
t

```

```
high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, colors.use = c("red", "grey"), do.hover =T, data.hover =c("ident", "res.1))
t

```

# Answer
> 2 votes
Since I don't think it is possible to have two color schemes in a single plot (one for expression and one for highlighting some cells) I will suggest other approach.

A possible option to do this could be to **change the shape** in the plot of the cells that you want to highlight and maintain the color code for expression.

I add a code that works for this purpose with the R object provided. Maybe there are easier ways to do that but this is working:

```
# Run PCA
seuset_h16 <- RunPCA(seuset_h16,
                     pc.genes = seuset_h16@var.genes)
# Run tSNE
seuset_h16 <- RunTSNE(object = seuset_h16,
                      dims.use = 1:10,
                      do.fast = TRUE)
# Find clusters of cells
seuset_h16 <- FindClusters(object = seuset_h16,
                           reduction.type = "pca",
                           dims.use = 1:10,
                           resolution = 0.6,
                           print.output = 0,
                           save.SNN = TRUE)
# Store output from FeaturePlot (do.return = T)
fp <- FeaturePlot(object = seuset_h16,
                  features.plot = c("DDB_G0273311"),
                  cols.use = c("grey", "orange"),
                  reduction.use = "tsne",
                  do.return = T)
# Add extra column in the data with the cells to be highligthed
fp$DDB_G0273311$data$cellToDisplay <- as.factor(ifelse(row.names(fp$DDB_G0273311$data) %in% c("s1.1", "s1.4", "s1.80"),
"group2",
"group1"))
# Modify color and shape in the plot
# lapply is needed since the object returned in fp is a list
lapply(fp, `+`, geom_point(aes(shape=cellToDisplay, colour = gene, size = 0.4)))

```

This is the output, with three cells highlighted as triangles and the rest as dots. You can play a little bit with shapes and sizes if you want a clearer display:

# Answer
> 0 votes
This updates the answer provided by Mack123456 to Seurat 3 and works in 3.1.2. There are some minor difference with colouring, but it is otherwise similar. The HoverLocator plot throws a warning, all about the plot colours and title, but the plot is fine and labels appear as expected.

```
load(file= "~/seuset_16.RData")

seuset_h16 <- UpdateSeuratObject(seuset_h16)

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(rownames(seuset_h16@meta.data), high)
Idents(seuset_h16, cells = high) <- "high"
Idents(seuset_h16, cells = low) <- "low"
seuset_h16$expr <- Idents(seuset_h16)
Idents(object = seuset_h16) <- "res.1"

t <- DimPlot(seuset_h16, pt.size = 2, shape.by = "expr")
t

t2 <- DimPlot(seuset_h16, pt.size = 2, cols = c("red", "orange","green","blue"))
t2 <- HoverLocator(plot = t2, information = FetchData(object = seuset_h16, vars = c("ident", "res.1")))
t2

```

---
Tags: r, single-cell, seurat, ggplot2
---


# Question
Title: How to Convert From Gene Location to SNP Location
I am new to this. I have an SNP array composed of digits \[0,1,2\] for chromosome 21. How do I find the start:end locations for the region that is "21q22.11" step by step?

# Answer
> 2 votes
You first need to get the coordinates of the "cytoband" 21q22.11, which I would search for using the UCSC table browser. You need the following settings: Clade: Mammal genome: Human assembly: group: Mapping and Sequencing track: Chromosome Band table: cytoBand region: position chr21 output format: BED

(you can leave other options as default)

That will give you a bed file of cytobands on chr21. You'll need to search (e.g. grep) for your cytoband: q22.11.

---
Tags: gene, snp
---


title: Wrong format .ped in PLINK

text:
# Question
Title: Wrong format .ped in PLINK
I know that I don't have the correct format in my .ped file for PLINK but I don't know how to convert it. My actual file is something like:

```
BSCAN BSCAN20016 0 0 2 0 0 1 1 1....
BSCAN BSCAN20017 0 0 2 0 0 1 1 1....
BSCAN BSCAN20018 0 0 2 0 0 1 1 1....

```

Where column 1 is FamilyID, column 2 is ID, column 3 is FatherID, column 4 MotherID, column 5 Sex and then the genotype of each individual. I only need to create binary files, so I tried to run:

```
plink-1.9-rc --file admixed --make-bed --no-pheno --cow --out admixed

```

but the program shows: Error: half-missing call in .ped file at variant 1, line 1.

Anyone know how can I change the format or recode it to make a correct input?

Thanks in advance!

# Answer
> 1 votes
So, columns 7 and 8 should be the genotype of variant 1 (each column is one of the alleles). You have them as 0 1, which means absent and major allele.

I think you're just missing a column. From the doc:

Column 6 is Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

Ideally you would know which SNPs are heterozygous based on your SNP calling and have such SNPs as 1 2, but if you have a homozygous organism or if you're working with data that's been flattened you can simply generate the .ped file with the same number twice for each SNP.

---
Tags: file-formats, plink
---


# Question
Title: How to use busco installed through conda?
I have installed busco on my Mac computer using `conda install --channel bioconda busco=3.0.2` into an environment called `busco3_0_2`. However, when I type `busco`, the command is not found.

Anyone knows what the name of the binary is? Usually when I install something with conda, the name of the command is the same as the name of the package. In this case, it does not seem to be so.

Note that the package installation finished successfully.

# Answer
From the README, it looks like the command-line entry point is not `busco` but `run_BUSCO.py`. This seems to agree with the bioconda recipe, where `run_BUSCO.py` and `run_busco` are listed among the commands.
> 4 votes

# Answer
The new version of busco (v4) has the user entry point called `busco` (it's still a python script). BUSCO v4 is available at conda, note it's a very fresh update.
> 1 votes

---
Tags: software-installation, conda, busco
---


# Question
Title: How I can push a function to act on my files
I have R code takes a vcf file returns a parsed .txt output.

I have 21 vcf files in a folder

I named my R code parse\_vcf\_alt1 function

# Answer
> 0 votes
By adding this to the script

```
args=(commandArgs(TRUE))

file_list <- list.files(path="/path/to/vcf/files/")
for (i in 1:length(file_list))
{
  x=file_list[i]
  fn=paste('test',x,sep = "")
  parse_vcf_alt1(x,fn)
}

```

But Also adding

```
df <- apply(vcf,2,as.character)
  fn=paste("C:/Users/44736/Videos/Captures/New folder/",fn,sep="")
  write.csv(df,fn)

```

to the main function I could get this

The problem is I don't know how to have a tab delimited file rather than csv also removing quote from the written output

# Answer
> 0 votes
Why don't you try using docopt which can make your R script "understand" command line parameters.

Doing that, you can create a simple bash loop to execute the script for all vcf files, such as:

`for files in "*.vcf" ; do script.R -i $files -o {output} ; done`

See this example here

Hope it helps, Felipe.

# Answer
> 0 votes
The upper part, where you try to use `lapply()` over your list of files with your custom function does not quite make sense to me: You declare your `files` but you do not use this in your `lapply()` call. I would expect to see something like \`lapply(files, your\_custom\_function).

Moreover, the error you get is results from `dplyr::rename()` not being able to find `#CHROM` because `readr::read_tsv()` could not read the header line properly. In my case the vcf was read as an n x 1 data frame! You could use `readr::read_delim()` to make sure you are reading the data right.

I would, however, would not try to read a `VCF` file as an ordinary data frame, there are packages for that, for example `[vcfR][1]`. Once you create your vcfR object, you can check the different slots:

> The vcfR object is an S4 class object with three slots containing the metadata, the fixed data and the genotype data.

---
Tags: r, loop
---


# Question
Title: Searching for proteins by primary and secondary structure
I need to find peptides/proteins with a certain sequence of amino acids in a certain type of secondary structure. For example, "ALA PHE GLY" in alpha-helix. Is there any tool that does that?

# Answer
Your request is rather specific, so there is no tool...

It is a straightforwards request though. I'd impose the requirement that the structures be solved and no predicted secondary structure. You can do it two ways if using Python.

## PDB

After downloading all PDB structures, in python with PyMOL (installable with conda, tricky without) you'd do something like

```
import pymol2, os
from multiprocessing import Pool

def find_motif(file):
    motif = 'AFG'
    try:
        with pymol2.PyMOL() as pymol
            pymol.cmd.load(file)
            seq = pymol.cmd.get_fastastr()
            if motif not in seq:
                return None
            d = {x: []}
            # str.find will cause issues if there is a gap also structure may not start at one.
            pymol.cmd.iterate('resn ALA+PHE+GLY', 'x.append((resi, resn, ss))', space=d)
            sequence = -1
            i = 0
            for residue in d['x']:
                if residue['resn'] = 'ALA':
                    sequence = 0 # lazy enum
                    i = int(residue['resi'])
                elif sequence != -1 :
                    if int(residue['resi']) == i+1 and residue['resn'] = 'PHE':
                        sequence = 1
                    elif int(residue['resi']) == i+2 and residue['resn'] = 'GLY':
                        return (file, i)
                    else:
                        sequence = -1
    finally:
        return (file, None)

if __name__ == '__main__':
    p = Pool(12)
    structure_files = [f for f in os.listdir('pdb') if os.path.splitext(f)[1] == '.pdb']
    verdicts = p.map(finder, structure_files)
    print([(f,v) for f,v in verdicts if v is not None])

```

The end result will be a list of PDB with where in the sequence (in that PDB code). (I wrote this in the reply box, so there may be a typo or two). PyMOL has a selection algebra called `pepseq`, which matches missing atoms, but I was a whimp and played it safe by being more verbose. Also, a for loop can be used instead of an iterate. But I think it is slower.

```
atoms = pymol.cmd.get_model(f"pepseq {target} and name CA")
for atom in atoms.atom:
      if atom.ss == ' ...

```

## Uniprot

In Uniprot, if a protein is crystallised, there will be the feature entries corresponding to the secondary structure. It is sheet or helix. No turns, pi-helices or 3\_10 helices. There is no Biopython method for parsing Uniprot XML. I wrote myself a parser and it is not too horrid thing to do —although I monkeypatched elementtree. Alternatively, `xmlschema.XMLSchema.to_dict` can be used to convert the XML to a dictionary. The advantage of Uniprot is that the residues ids will not require correcting. The disadvantage is parsing XML.

My parser is specific for my application and it is poorly documented (so would not recommend it) but but... I can share it and human data. I'd be used like:

```
 import os
 from michelanglo_protein import ProteinCore
 ProteinCore.settings.startup(data_folder='protein-data')
 p = ProteinCore().gload(file=os.path.join('gpickle','Q86V25.gpz'))
 #print(p.asdict())
 if 'helix' in p.features:
        n = p.sequence.find( ...
        for ... p.features['helix'] #('helix', 'turn', 'strand')
        if x < n < y ...

```
> 3 votes

---
Tags: protein-structure, secondary-structure
---


# Question
Title: Downloading dataset from SRA (SOLiD Platform)
Trying to download colorspace data from SRA, but getting an error

`abi-dump -A SRR1657115.sra`

abi-dump.2.9.6 err: item not found while constructing within virtual database module - the path 'SRR1657115.sra' cannot be opened as database or table.

# Answer
You have to use the accession name without the `.sra` suffix.

`abi-dump -A SRR1657115`
> 1 votes

---
Tags: ngs
---


title: Removing rows containing some strings

text:
# Question
Title: Removing rows containing some strings
I have a big data frame

```
> dim(anno_maf)
[1] 449892    132
>

```

In chromosome column I have chr 1 to 22 and another contigs like

```
> unique(anno_maf$Chromosome)
 [1] "chr1"                  "chr6"
 [3] "chr16"                 "chr17"
 [5] "chr20"                 "chr2"
 [7] "chr3"                  "chr4"
 [9] "chr14"                 "chr19"
[11] "chr5"                  "chr10"
[13] "chr9"                  "chr12"
[15] "chr13"                 "chr11"
[17] "chr22"                 "chr7"
[19] "chr15"                 "chr8"
[21] "chr18"                 "chr21"
[23] "chrX"                  "chrY"
[25] "chr4_gl000194_random"  "chr17_gl000205_random"
[27] "chrUn_gl000241"        "hs37d5"
[29] "chrUn_gl000219"        "chrUn_gl000234"
[31] "chr1_gl000191_random"  "chrUn_gl000211"
[33] "chrUn_gl000224"        "chrUn_gl000225"
[35] "chr17_gl000203_random" "chrUn_gl000212"
[37] "chrUn_gl000243"        "chrUn_gl000214"
[39] "chrM"                  "chr1_gl000192_random"
[41] "chr7_gl000195_random"  "chrUn_gl000232"
[43] "chr4_gl000193_random"  "chr19_gl000208_random"
[45] "chrUn_gl000226"        "chrUn_gl000218"
[47] "chr9_gl000199_random"  "chrUn_gl000217"
[49] "chrUn_gl000229"        "chrUn_gl000216"
[51] "chrUn_gl000231"        "chr9_gl000198_random"
[53] "chr17_gl000204_random" "chrUn_gl000220"
[55] "chrUn_gl000235"        "chr11_gl000202_random"
[57] "chrUn_gl000222"        "chrUn_gl000240"
[59] "chrUn_gl000233"        "chrUn_gl000230"
[61] "chrUn_gl000213"        "chrUn_gl000238"
[63] "chr19_gl000209_random" "chrUn_gl000237"
[65] "#CHROM"

  How I can only keep rows having chr1 to chr22 and removing the rest?

```

# Answer
You can create a vector consisting of values that you would like to keep (or discard):

```
to_keep <- paste0("chr", seq(1,22))

> to_keep
 [1] "chr1"  "chr2"  "chr3"  "chr4"  "chr5"  "chr6"  "chr7"  "chr8"  "chr9"  "chr10" "chr11"
[12] "chr12" "chr13" "chr14" "chr15" "chr16" "chr17" "chr18" "chr19" "chr20" "chr21" "chr22"

your_df[your_df$Chromosome %in% to_keep,]

```

The `data.table` way, which should be faster for large data, would be:

```
your_df[Chromosome %in% to_keep]

```

I don't know what your goal is but you might want to keep chromosomes X and Y as well.
> 1 votes

---
Tags: r
---


# Question
Title: How do I find identical sequences in a FASTA file?
I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in Biopython for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).

Using the following example input

```
>prot_id_1
MAWIGLDISKLFVENSRDAAA
>prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_5
MAWIGLDISKLFVENSRDAAA
>prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_7
MAWIGLDISKLFVENSRDAAA

```

I currently do it as follows:

```
from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) > 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')

```

which gives me the desired output:

```
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS

```

Is there a smarter way of accomplishing this?

# Answer
> 8 votes
The trick would be to swap the key in the dictionary to be the sequence itself.

Also I would recommend using a different separator that "\_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.

```
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")

```

# Answer
> 4 votes
Just as a small variation to @Bioathlete's answer in case you want to write the fasta using Biopython (e.g. to add names and description):

```
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)

# this creates a generator; if you need a physical list, replace the outer "(", ")" by "[" and "]", respectively
final_seq = (SeqRecord(Seq(seqi, IUPAC.protein), id="|".join(gi), name='', description='') for seqi, gi in dedup_records.items())

# write file
SeqIO.write(final_seq, 'test_output_final.fasta', 'fasta')

```

This will give the desired output:

```
>prot_id_1|prot_id_5|prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_2|prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

```

# Answer
> 2 votes
using **datamash** ans some awk

```
$ cat input.fa | paste - -  |\
cut -c2- | sort -t $'\t' -k2,2 |\
datamash  collapse 1 -g 2  | tr "," "_" |\
awk '{printf(">%s\n%s\n",$2,$1);}'

>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

```

# Answer
> 1 votes
This can also be done quite easily with awk:

`awk '{if(NR%2){printf("%s\t",$1)}else{print}}' filename | sort -k2,2| awk '{if($2==SEQ){gsub(">","",$1);ID=ID"|"$1} else{if(SEQ!=""){printf("%s\n%s\n", ID,SEQ);}SEQ=$2;ID=$1;} }END{printf("%s\n%s\n", ID,SEQ)}'`

The first `awk` converts the fasta file to a tab separated file with format `ID\tSequence`, which is then sorted by sequence by `sort`. The last `awk` goes through the sorted file looking at the sequences: if the sequence in the current line is the same as that in the previous line, it appends its id to the ID variable. Otherwise (i.e. we have encountered a new sequence), it prints the previous concatenated ids and the corresponding sequence.

# Answer
> 1 votes
Here is my free program on Github **Sequence database curator** (https://github.com/Eslam-Samir-Ragab/Sequence-database-curator)

It is a very fast program and it can deal with:

1. Nucleotide sequences
2. Protein sequences

It can work under Operating systems:

1. Windows
2. Mac
3. Linux

It also works for:

1. Fasta format
2. Fastq format

Best Regards

# Answer
> 1 votes
If you have a very big fasta file(60 GB) why use SeqIO.parse and write the output file using biopython?

Why not use SimpleFastaParser instead (you get a tuple(name, seq)) and you dont need to wait so long as using the parser and using Seqrecords.

What do you think about it guys?

Thanks

Paulo

Kind of:

```
from Bio.SeqIO.FastaIO import SimpleFastaParser
from collections import defaultdict

dedup_records = defaultdict(list)
for (name, seq) in SimpleFastaParser(filename):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[seq)].append(name)
with open(filename, 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")

```

---
Tags: fasta, python, filtering
---


title: Merging two data frames

text:
# Question
Title: Merging two data frames
I have a data frame like with more samples like

```
> head(a[1:2,1:4])
             ID               ids pre.or.post.treatment
1 OCCAMS/AH/120 LP6005409-DNA_F01                   pre
2 OCCAMS/AH/126 LP2000325-DNA_A01                   pre
            TRG
1 non-responder
2 non-responder
> dim(a)
[1] 160   4
>

```

Another data frame with part of my samples but more characteristics for each sample like

```
> dim(b)
[1] 60 50

> head(b[1:2,1:6])
             ID       TRG pre.or.post.treatment          ID.1
1 OCCAMS/AH/212 responder                   pre OCCAMS/AH/212
2 OCCAMS/AH/279 responder                   pre OCCAMS/AH/279
  DiagnosisDate FE.LastSeenDate
1      8/8/2013        6/6/2014
2     6/10/2014       11/8/2018
>

```

How I can merge these data frame getting a new data frame with my all 160 samples columns of both a and b

# Answer
> 2 votes
Assuming your column ID in both dataframes `a` and `b` shared same IDs, and as your dataframe are not particularly heavy, you can have the use of the function `left_join` from `dplyr` package as this:

```
library(dplyr)
DF <- left_join(a,b, by = "ID")

```

# Answer
> 3 votes
@dc37's dplyr solution is quite elegant, but you can also use base R's merge:

```
merge(x = df_a, y = df_b, by = "ID", all.x = TRUE, all.y = TRUE)

```

You can tweak various parameters, such as use `by.x` and `by.y` instead of just `by`, change `all.x` and `all.y` based on requirements to get your results.

This is a pretty common R question by the way, and if you had googled "merging data frames", you'd see a lot of options.

---
Tags: r
---


# Question
Title: how can i calculate p value of miRNA raw reads counts using R?
working on lung cancer project and need to find out which genes are deferentially expressed and which are not. for that i a have downloaded miRNA data set from GEO id is GSE110907. And the file i got has to variable gene\_id and read counts. from that i have to derive p value from that data/file.

# Answer
I would start by re-formulating the question in this post: Your end-goal seems to the differential expression analysis (DEA) of your samples (and you have to have to levels for this like treatment-control, mutant-WT, etc ...). P-values are context-dependent and should be calculated based on your research question/hypothesis.

DEA has been addressed by many researchers and `R` has a lot of packages to help you with this; `DESeq2`, `edgeR`, `limma`, etc ...

I would go to the `Biconductor` website, search for one of the three above, and go through their rich **vignettes** and would come back here for **specific** questions.
> 1 votes

---
Tags: r, p-values
---


title: Chaning string in rows

text:
# Question
Title: Chaning string in rows
I have a too big data frame

```
> dim(anno_maf)
[1] 2635299     132
>

> head(anno_maf[1:3,1:4])
                     Tumor_Sample_Barcode Chromosome Start_Position End_Position
1: LP2000107-DNA_A01_vs_LP2000102-DNA_A01       chr5      137536800    137536803
2: LP2000330-DNA_A01_vs_LP2000316-DNA_A01      chr11       64012229     64012231
3: LP2000330-DNA_A01_vs_LP2000316-DNA_A01      chr17       41720917     41720917
>

> unique(anno_maf$Tumor_Sample_Barcode)
  [1] "LP2000107-DNA_A01_vs_LP2000102-DNA_A01" "LP2000330-DNA_A01_vs_LP2000316-DNA_A01" "LP2000332-DNA_A01_vs_LP2000323-DNA_A01"
  [4] "LP2000333-DNA_A01_vs_LP2000324-DNA_A01" "LP6005334-DNA_A02_vs_LP6005333-DNA_A02" "LP6005334-DNA_C01_vs_LP6005333-DNA_C01"

```

I want to remove everything thing from `vs` onward in `anno_maf$Tumor_Sample_Barcode`, for example leaving `LP2000107-DNA_A01` from `LP2000107-DNA_A01_vs_LP2000102-DNA_A01` but I don't know how to do that

Any help?

# Answer
> 3 votes
`sub()` is the way to go.

Your pattern is your regular expression: Check for "\_vs\_" followed by any character (provided with `.`) for an undefined length (provided with `*`).

```
my_text <- "LP2000107-DNA_A01_vs_LP2000102-DNA_A01"

sub(pattern = "_vs_.*",
    replacement = "",
    x = my_text)

[1] "LP2000107-DNA_A01"

```

---
Tags: r
---


title: Normalization using parallel R script

text:
# Question
Title: Normalization using parallel R script
I need to implement a parallel version of an R script because on one core it takes forever to execute. The function that takes the most of time is `scran::quickCluster()` applied to an integer sparse matrix (`UMI_count`) that is `26000x66000`. The code is presented below:

```
library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, cluster=clusters)

```

I did not find much on parallel versions of `quickCluster()` function, so will it actually execute in parallel if I set the number of cores to more than one? I am trying to do normalization of heterogeneous dataset with multiple cell types, that is why I am using `quickCluster()`.

And I just need to speed it up to be able to compute it.

# Answer
> 2 votes
The `scran::quickCluster` method has a `BPPARAM` (i.e., `BiocParallelParam`), so if one provides something like

```
clusters <- quickCluster(sce, BPPARAM=BiocParallel::MulticoreParam())

```

it will use all available cores (a specific number of cores can also be given). See the intro doc to `BiocParallel` for info on other options.

---
Tags: r, scran
---


# Question
Title: Stacked barplot for single cell analysis
I want to create a stacked bar-graph with different cell cycles for each cell type within each condition. I have uploaded the file for it. Don't know how to go about it as now I have another condition (cell cycle). I was able to do the bar graph comparing cell cycle per condition but now I want to compare cell type and cell cycle in each condition

output of dput(head(all.combined@metadata))

```
structure(list(orig.ident = c("treated", "treated", "treated",
    "treated", "treated", "treated"), nCount_RNA = c(1892, 307, 1348,
    3699, 4205, 4468), nFeature_RNA = c(960L, 243L, 765L, 1612L,
    1341L, 1644L), percent.mt = c(0.211416490486258, 1.62866449511401,
    4.45103857566766, 4.4065963773993, 0.0713436385255648, 3.87197851387645
    ), RNA_snn_res.0.5 = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
    "13", "14", "15", "16", "17", "18", "19"), class = "factor"), seurat_clusters = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"), S.Score = c(0.476893835992198, -0.0200784617568548, -0.0335915198305002, -0.0247184276246385, 0.010785196602457, 0.0190008903712199), G2M.Score = c(0.204441469200986, 0.173804859670862, -0.0313235510969097, -0.0376796363661889, -0.0559526905696905, -0.0122031631356698), Phase = structure(c(3L, 2L, 1L, 1L, 3L, 3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(7L,7L, 1L, 4L, 7L, 9L), .Label = c("Fibroblast", "T cell", "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("treated_AAACGCTAGCGGGTTA-1", "treated_AAAGGTAAGTACAGAT-1", "treated_AAAGTGAGTTTGATCG-1", "treated_AAATGGACAAAGTGTA-1",
    "treated_AACAAAGGTCGACTTA-1", "treated_AACAGGGTCCTAGCCT-1"), class = "data.frame")

```

output of dput(tail(all.combined@metadata))

```
structure(list(orig.ident = c("untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), nCount_RNA = c(901, 823,
1184, 1835, 1147, 1407), nFeature_RNA = c(482L, 479L, 649L, 1043L,
604L, 709L), percent.mt = c(1.77580466148724, 2.91616038882138,
4.22297297297297, 3.86920980926431, 2.0052310374891, 4.05117270788913
), RNA_snn_res.0.5 = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
"13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    seurat_clusters = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
    "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",
    "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    S.Score = c(-0.0320858200243315, 0.0304725660342869, 0.0215996091745327,
    0.0384166213301423, 0.144956251122548, -0.0242770509986111
    ), G2M.Score = c(0.0904224391544142, 0.050148242050667, -0.0178041670730754,
    -0.0112596867977946, -0.0519554524339088, -0.0136533184257381
    ), Phase = structure(c(2L, 2L, 3L, 3L, 3L, 1L), .Label = c("G1",
    "G2M", "S"), class = "factor"), old.ident = structure(c(5L,
    5L, 5L, 5L, 5L, 5L), .Label = c("Fibroblast", "T cell", "Macrophage",
    "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
    "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("untreated_TTTGGTTGTCTAATCG-18",
"untreated_TTTGGTTTCCCGAGGT-18", "untreated_TTTGTTGAGAACTGAT-18",
"untreated_TTTGTTGAGCTCGGCT-18", "untreated_TTTGTTGAGTGCCTCG-18",
"untreated_TTTGTTGCACGGTGCT-18"), class = "data.frame")

```

Been using this code to generate the previous graph.

```
ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

```

# Answer
Based on your example, you can use count for each phase of cell cycle for each cell types in each condition. Here, I'm using `dplyr` package to do that but most likely, you could have the same output using various methods:

```
DF <- rbind(dfhead, dftail)

library(dplyr)
DF_Count <- DF %>%group_by(orig.ident,Phase,old.ident) %>%
  count() %>%
  ungroup() %>%
  group_by(orig.ident,old.ident) %>%
  mutate(Freq = n/sum(n)*100)

# A tibble: 8 x 5
# Groups:   orig.ident, old.ident [5]
  orig.ident Phase old.ident      n  Freq
  <chr>      <fct> <fct>      <int> <dbl>
1 treated    G1    Fibroblast     1 100
2 treated    G1    Stellate       1 100
3 treated    G2M   Tumor          1  33.3
4 treated    S     Tumor          2  66.7
5 treated    S     Mast cell      1 100
6 untreated  G1    Acinar         1  16.7
7 untreated  G2M   Acinar         2  33.3
8 untreated  S     Acinar         3  50

```

As you can see, `DF_Count` has the frequency for each cell types of each phase of the cell cycle in function of the condition. We can use `DF_Count` to get the following plot. Using `facet_wrap`, you can create two panels based on the condition column and thus represent the cell cycle of each cell types in function of the treatment condition:

```
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

```

Does it look what you are expecting ?

NB: `DF` here is the addition of your head and tail of your `all.combined` dataframe.
> 1 votes

---
Tags: single-cell, ggplot2
---


title: Moving file based on their names

text:
# Question
Title: Moving file based on their names
I have a list of vcf files; I also have a list of names in a txt file like

```
LP6005409-DNA_F01
LP2000325-DNA_A01
LP6005409-DNA_E02
LP6005500-DNA_C03

```

What I have in this txt files are part of name of each vcf files

I want to move .vcf files having these in their name to a separate folder

I done

```
[fi1d18@cyan01 TRG]$ for i in *.vcf
> do if grep -q $i 1.txt;
> cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
> done
-bash: syntax error near unexpected token `done'
[fi1d18@cyan01 TRG]$

```

How I can complete this task?

Thank you

# Answer
First of all you are getting an error because you are missing a `then` in the third line of your code and then a `fi` to close the `if`. It should be:

```
for i in *.vcf
do if grep -q $i 1.txt; then
    cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
fi
done

```

And I guess the one-liner below achieves what you are trying to do:

```
while IFS= read -r line; do mv "$line"*.vcf some-directory; done < file-list.txt

```
> 4 votes

# Answer
`for i in $(cat text.txt); do cp -v *${i}*vcf mydir; done`
> 1 votes

---
Tags: bash, linux
---


# Question
Title: Is it possible to obtain RefSeq transcripts from Ensembl database for GRCh37 to match Ensembl Genes?
I know that with GRCh38, it is possible to extract the matching RefSeq transcripts for Ensembl genes/transcripts/proteins using Biomart from Ensembl.

However, for GRCh37, this option does not exists. Using the following MySQL command returns Ensembl transcript IDs for matching RefSeq transcripts:

```
mysql -u anonymous -h ensembldb.ensembl.org -A -P 3337 -D homo_sapiens_core_99_37 -e "SELECT transcript.stable_id, xref.display_label FROM transcript, object_xref, xref,external_db WHERE transcript.transcript_id = object_xref.ensembl_id AND object_xref.ensembl_object_type = 'Transcript' AND object_xref.xref_id = xref.xref_id AND xref.external_db_id = external_db.external_db_id AND external_db.db_name = 'RefSeq_mRNA';"

```

Also, I was trying to play around with the databases schema but I am not able to find any solution to problem.

Ultimately, I wish to obtain a table like the following:

```
Gene stable ID      Transcript stable ID    Protein stable ID   RefSeq match transcript
ENSG00000120329     ENST00000239451         ENSP00000239451     NM_031947.4
ENSG00000164219     ENST00000419445         ENSP00000404676     NM_005023.4
ENSG00000150594     ENST00000280155         ENSP00000280155     NM_000681.4
ENSG00000113578     ENST00000337706         ENSP00000338548     NM_000800.5

```

Thank you in advance.

# Answer
The RefSeq match option in BioMart is from the Matched Annotation from NCBI and EBI (MANE) collaboration between RefSeq and Ensembl. It has only been calculated for the up-to-date gene annotation on GRCh38 so cannot be obtained on GRCh37. You can get mapping from Ensembl to RefSeq transcripts through BioMart as RefSeq mRNA ID (refseq\_mrna in R) but this is not a perfect match like the MANE, it is a mapping based on sequence similarity and similar genomic location, and there can be mismatches between them.
> 4 votes

# Answer
Try the GRCh37 version of BioMart, and the External reference "RefSeq mRNA ID" / "RefSeq ncRNA ID":

https://grch37.ensembl.org/biomart

From the Perl code provided by BioMart, it looks like the field name should be "refseq\_mrna" or "refseq\_ncrna". I'm not sure how that translates to the SQL code.

I feel compelled to add a caution to this, though. GRCh37 is a very old build, having been *replaced* by GRCh38 in December 2013. If possible, please find a way to use GRCh38 for your analysis instead, and don't create any *new* software that depends on GRCh37. For example, gnomAD v3 supports GRCh38 (see bioconductor, or here).
> 1 votes

---
Tags: ensembl, refseq
---


# Question
Title: Is there a point in recalibration of scores for variant calling?
The most variant calling pipeline GATK include a Base Quality Score Recalibration (BQSR) which requires a list of known variants. Recently, some work has been done for reference-free recalibration of scores as well: Lacer and atlas, which is motivated by making the most for aDNA and low coverage datasets.

The importance for aDNA is explained in this lecture, but it is not clear to me if / how is important BQSR is for fresh DNA samples with decent (\>15x) coverage. Especially when I work with non-model organisms and I can not simply use the standard tools.

How big an impact does recalibration of scores have on variant calling? Is there a rule of thumb for which it is / it is not worth the effort?

# Answer
I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few points at most, but it's certainly noticeable.

GATK recommends BQSR for both genome and exome data, which is normally much higher than 15x.
> 5 votes

# Answer
That's a good question.

I'd say that you don't need to bother with variant recalibration for

* low number of samples (e.g., just two trios); I could not get GTAK recalibration of variant scores to work anyway
* high-coverage samples (e.g., X Ten genomes with 30x coverage) where the DNA samples themselves are of high, comparable quality and have been sequenced with consistent technology.

Generally, it is my impression that a lot of the thoughts and advanced statistical models built into GATK come from the earlier phases of the 1000 Genomes project. This means (1) low-coverage, (2) different coverage genomes (3) sequenced with varying technology versions by (4) different samples and (5) population sequencing.

If you are in a clinical setting where you do 30x sequencing on X Ten platforms only anyway, then variant recalibration will probably not help you that much.

On the other hand, if you are integrating many data sets from different data centers and machine versions etc., variant recalibration might be worth a shot.

A good check would be looking at genotype quality distributions and other variant/quality related metrics before and after recalibration.

Anyone: please correct me if I'm wrong!
> 2 votes

# Answer
Ideally these BQSR methods were made keeping in mind how technical errors will actually screw up the base quality calls and when the machines were still more on development phase while being used for the 1000G project. As of now machines are more powerful and strong where it will be unlikely to use it but still we use with listed SNPs to find the covariates and build a model around the data using the information with machine learning tricks to improve the quality of those base calls. Ideally it should be more appropriate when old machines from Illumina or other standard companies are being used but with new machines which are much powerful and having high throughput they should tend to go down. I do not recall if such tests have been made but obviously I know new sequencing machine always make such tests to show that they have reduced such errors but still recommend such BQSR for variant calls. Now the problem is the list of SNPs, this to me is the real problem since the list we use is far from being Gold standard and if that is not properly taken care of everything we infer about quality is still shaky. This link is pretty informative but it's an old one. I would really see improvements with new sequencers. However very less people care about such tests in academic research and also translational lab will really not invest time and money on such unless the facility has some bioinformaticians who always does such testing while buying a new sequencer for the institute. In terms of clinical genomics for finding variants I reckon most powerful and up-to-date sequencers should be used but not sure if they still use BQSR and if so what is the list they use to build model of covariation around the data.
> 1 votes

# Answer
In case BQSR is not an option (i.e. non-model organisms) it would be best to use some internal control sequence such as PhiX for illumina platform. Although this is supposed to be common practice some facilities ignore it. In principle, the machines should use these sequences as reference so that the scoring would be more accurate. In my experience the first 10-15 bases of the illumina reads always had lower quality. This can easily be seen in the nucleotide distributions. I would advise trimming the first 10-15 bases and quality based end-trimmingiIf the quality of the individual reads are important, such as low coverage resequencing or de-novo genome assembly applications.
> 0 votes

---
Tags: ngs, variant-calling
---


# Question
Title: Making a boxplot or violinplot for several dataframe with different number of samples
I have 3 data frames for treatment of 3 different drugs and in each of them I have the number of mutation types like insertion, deletion, SNP and total of mutations for each patients. In each group I have different number of patients like below

```
> head(dat1)
                patient  DEL  INS   SNP MNP total
1:    LP6008337-DNA_H06  927  773 40756   0 42456
2:    LP6008334-DNA_D02 1049  799 31009   0 32857
> dim(dat1)
[1] 21  6
>

 head(dat2)
                Patient   DEL  INS   SNP MNP total
1:    LP6008031-DNA_E01 13552 3374 62105   0 79031
2:    LP6005500-DNA_G01   539  500 43451   0 44490
> dim(dat2)
[1] 33  6
>

> head(dat3)
                Patient   DEL   INS    SNP MNP  total
1:    LP6005935-DNA_F03 39168 16739  58095   0 114002
2:    LP6008269-DNA_D08   849   910 103501  11 105271

> dim(dat3)
[1] 106   6

```

I want to show if the number of each mutational category and total is different for each drug by boxplot, violinplot or any type of visualization

For example this picture is very nice

Can you help me?

# Answer
To help us help you, please make your data available with `dput()` next time.

To achieve your goal with `ggplot2`, you would need all of your data in one data frame and in the "long format".

```
dat1 <- "
Patient DEL INS SNP MNP total
LP6008337-DNA_H06 927 773 40756 0 42456
LP6008334-DNA_D02 1049 799 31009 0 32857
"

dat2 <- "
Patient DEL INS SNP MNP total
LP6008031-DNA_E01 13552 3374 62105 0 79031
LP6005500-DNA_G01 539 500 43451 0 44490
"

dat3 <- "
Patient DEL INS SNP MNP total
LP6005935-DNA_F03 39168 16739 58095 0 114002
LP6008269-DNA_D08 849 910 103501 11 105271
"

# read data and add a column indicating data frame of origin (drug in this case)
dat1 <- readr::read_delim(dat1, delim = " ")
dat1$drug <- "drug1"
dat2 <- readr::read_delim(dat2, delim = " ")
dat2$drug <- "drug2"
dat3 <- readr::read_delim(dat3, delim = " ")
dat3$drug <- "drug3"

# bind rows
dat <- rbind(dat1, dat2, dat3)

> dat
# A tibble: 6 x 7
  Patient             DEL   INS    SNP   MNP  total drug
  <chr>             <dbl> <dbl>  <dbl> <dbl>  <dbl> <chr>
1 LP6008337-DNA_H06   927   773  40756     0  42456 drug1
2 LP6008334-DNA_D02  1049   799  31009     0  32857 drug1
3 LP6008031-DNA_E01 13552  3374  62105     0  79031 drug2
4 LP6005500-DNA_G01   539   500  43451     0  44490 drug2
5 LP6005935-DNA_F03 39168 16739  58095     0 114002 drug3
6 LP6008269-DNA_D08   849   910 103501    11 105271 drug3

# convert data to long format, ggplot() expects that
dat_long <- tidyr::pivot_longer(dat, cols = 2:6)

> dat_long
# A tibble: 30 x 4
   Patient           drug  name  value
   <chr>             <chr> <chr> <dbl>
 1 LP6008337-DNA_H06 drug1 DEL     927
 2 LP6008337-DNA_H06 drug1 INS     773
 3 LP6008337-DNA_H06 drug1 SNP   40756
 4 LP6008337-DNA_H06 drug1 MNP       0
 5 LP6008337-DNA_H06 drug1 total 42456
 6 LP6008334-DNA_D02 drug1 DEL    1049

ggplot(dat_long,
                aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug)

```

EDIT: setting `scales` argument of `facet_wrap()` to "free" could be helpful

```
ggplot(dat_long,
       aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug, scales = "free")

```
> 2 votes

# Answer
If you use a violin plot instead and ggplot2 to make it then:

```
geom_violin(...not relevant stuff..., scale="count")

```

will scale each violin plot according to the number of actual entries in it. I often find this nice so small groups look super skinny and people don't focus on random changes in them so much.
> 1 votes

---
Tags: r, visualization, ggplot2, wgs
---


title: RIP-seq analysis?

text:
# Question
Title: RIP-seq analysis?
Given an experiment consisting of an input (baseline RNA) and IP (pulldown to find RNAs bound to certain protein of interest)... Is a DE analysis performed over the RNA-seq data from the samples (lets say with EdgeR or DESEQ2) suitable to reveal the preferentially bound RNAs? What other software tools would you recommend?

# Answer
You need to be careful of terminology. To me, a RIP-seq experiment involves a pull down, followed by a RNAseq library prep. Thus, the whole transcript is captured, not just the binding site of the protein (as in CLIP-Seq, HITS-CLIP, PAR-CLIP, iCLIP or eCLIP). Thus "peak-callers" whether they be designed for calling protein-DNA or protein-RNA binding sites are not suitable as the signal will not be punctuate or peaky. Of the methods mentioned by @Devon Ryan, only RIPSeeker seems setup to deal with this sort of data.

In some of the papers that talk about calling RNA-protein interactions they lump RIP-seq in with the CLIP techniques (see piranha and ASPeak papers referenced by @Devon Ryan). They appear to be talking about a technique where the RNA is fragmented before it is pulled down. Thus you would capture peaky binding sites (actually, the protocol would be remarkably similar to CLIP).

You should work out which of these you have.

As for using DESeq/EdgeR etc, in the RIPSeeker paper they say:

> Furthermore, programs for de novo transcript assembly followed by differential expression (DE) analysis, such as the Cufflinks/Cuffdiff suite (15,16), and for DE on a set of known transcripts, such as DESeq (17), may appear applicable to RIP-seq analysis. Unlike peak-calling strategy, however, the transcript-based methods assume the full transcriptome being sequenced at a fairly deep coverage (as usually the case in RNA-seq) and thus may be sensitive to background noise typical to the IP-based protocols, which is due to both the non-specific RNA interactions with a protein of interest and the non-specific RNA input from the pull-down of the (mutant) control (Supplementary Figures S3 and S4).

I don't know to what extent RIPSeeker performs better than the naive DESeq approach as its not included as one of the comparators in the paper.
> 4 votes

# Answer
You can try doing standard differential expression, but I worry that the between-sample normalization will work poorly. Personally, I would do peak calling instead, followed by diffBind. You have a few tools to choose from when it comes to this. In the past, I've rolled my own methods for this using MACS2 and genomic alignments (I then converted those to bedGraph files where entries are transcripts, so peak calling is then in transcript-coordinates).

It's probably easier, though, to use premade software. There are a number of packages out there for RIP-seq. These include piranha, RIPSeeker and ASPeak. Particularly if you're just starting out, you'd be well served with using one of those tools.
> 5 votes

# Answer
I've got the RIP-seq analysis outputs. We used an isotype control antibody VS antibody against an RBP to pull down corresponding RNAs!

I did the following pipeline: Hisat2 (mapping the reads after QC) -\> htseq-count (getting read counts) -\> DESeq2 (getting enriched RNAs in test VS isotype control) -\> volcano plot shows up-regulated and down-regulated genes in the test sample.

Logically, I should have seen just enriched RNAs (up-regulated genes), but surprisingly there are some RNAs down-regulated in the test sample (or up-regulated in the isotype control sample)!!!!! How can we explain and interpret this???? it seems some genes are pull down in isotype control sample as well. Are these background? Is my computational pipeline right?

BTW, it's worth to mention that all enriched RNA (called up-regulated genes in DESeq2 pipeline) are significantly relevant and meaningful to my experiment. However, I am not sure how much confident I can be on this output? And how should I explain down-regulated RNAs? because normally we just expect enriched RNAs in RIP-seq experiment.

Any thoughts?
> 0 votes

---
Tags: rna-seq, software-recommendation
---


# Question
Title: How do I get the compartment Id and volume from SBML using libSBML Python?
I searched the Python documentation for libSBML but couldn't find a staightforward answer to this question. For a given SBML model, how can I use Python libSBML to retrieve the values of the `id` and `volume` attributes for a given compartment?

# Answer
> 1 votes
There aren't any examples in the python docs on how to do this. The best bet is to search the C/C++ docs which are more extensive. There is a python example in the docs which shows how to retrieve the number of compartments in a model but it doesn't go further. The code below shows how you can retrieve the information requested, assuming you have a model called testmodel.xml:

```
import libsbml

document = libsbml.readSBML("testmodel.xml")

if document.getNumErrors() > 0:
   print("Encountered the following SBML errors:" + "\n")
   document.printErrors()
else:
   model = document.getModel()

   nCompartments = model.getNumCompartments()

   print ("Number of compartments: ", nCompartments)
   for i in range (nCompartments):
       comp = model.getCompartment(i)

       print ("Id = ", comp.getId())
       print ("Volume = ", comp.getVolume())

```

---
Tags: sbml
---


# Question
Title: Produce a single sequential FASTA sequence out of BAM
I'm having problems properly looking for a solution because I'm a layman in Bioinformatics not familiar with the terminology. I'm hoping you can nudge me in the right direction, please! Thank you very much.

# What I want in the end

A FASTA file containing a single sequential sequence (chromosome 1-22..) of the human genome.

# What I have

* 23gb of "clean\_data" which are .fq.gz files.
* 21gb of "result\_alignment" which is .bam+.bai
* "result\_variation" indel.vcf.gz and snp.vcf.gz

# What I tried

I used `samtools fasta reads.bam > reads.fasta` which results in many lines like so:

```
>V100005418L4C001R001322578/1
AAATAGCCCACAC...
>V100005418L4C001R006084507/1
GCCATAAAGCCTA...
...

```

But these include many many duplicate lines and the reference genome from the human genome project GRCh38 is only 3GB in size, so I thought that this FASTA sequence needs to be properly "aligned" somehow?

Am I thinking correctly here?

I'm really confused with all of the terminology and don't know exactly what I need to look for. Could someone kindly put the process into simpler words for me?

# Answer
What I wanted to do is called "Consensus Sequence". Two steps were needed:

```
$ bcftools mpileup -Ou -f ref.fa input.bam | bcftools call -Ou -mv | bcftools norm -f ref.fa Oz -o output.vcf.gz

```

After that you can create the consensus:

```
$ bcftools consensus -f ref.fa output.vcf.gz > out.fa

```
> 2 votes

# Answer
If you just want the human chromosomes in a FASTA format, why not downloading directly the individual chromosome files (`chr*.fa.gz`)?

If you don't rely on the "official" releases (for whatever reason), then what you need is first to assemble your reads into genomes, possibly (but not surely) going down to a chromosome-level resolution. Different approaches for assembly are possible/available, mainly depending on (1) what do you biologically need to answer and (2) the type of sequencing data you have - we'll need more details to be able to help you in this direction.

The `samtools fasta` utility will "just" convert your reads (likely given from a sequencing facility) from an alignment (BAM) format into a reads (FASTA) format, there is no concept of chromosomes here. You won't get assembled chromosomes this way, just extract the sequence in FASTA format from your alignment file.
> 3 votes

---
Tags: sam, sequence-alignment, ngs, phylogenetics, genome-sequencing
---


# Question
Title: What is the recognized way to reference accession data (.sra files) from NCBI as a URI?
I would like to assign a URL/URI to accession runs (i.e. SRR IDs) from the NCBI short read archive.

In program listings, I used to list them as `ftp://` URLs, because they were available on their servers as such, e.g. `SRR3418242`:

```
ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR341/SRR3418242/SRR3418242.sra

```

But this approach is deprecated, plus, not all runs are available this way anymore. The now recommended approach is to use sra-tools to download the data. This tool takes the SRR Ids directly, i.e. `prefetch "SRR3418242"`, and resolves them to one of many possible locations (home directory, http, etc).

Is there an official URI scheme that can be used to reference accession runs from the short-read-archive ? I would like to disambiguate them from other types of files and URLs (S3 keys, http, file). For file manifests, spreadsheets, and publications, is there something similar to `srr://SRR3418242`?

# Answer
> 1 votes
The nomenclature is no longer being followed as NCBI is moving all the files to S3/Google cloud. pysradb allows you to fetch the metadata and get the URLs for each individual SRR. See examples in this notebook.

You can also download an entire project by doing:

```
pysradb download -p <SRP>

```

---
Tags: phylogenetics, sratoolkit
---