text:
RATT works on example bacterial sequence and other bacterial genome but not on C. elegans genome and annotation
I am trying to use a tool called RATT described in this paper and available here. I ran the example using the following command:
$RATT_HOME/start.ratt.sh embl F11.fasta F11 Strain
embl is a folder containing the file Tb_H37Rv.embl, which is an annotation file in EMBL format of the Mycobacterium tuberculosis genome. F11.fasta contains the genome in FASTA sequence of the strain F11 of Mycobacterium tuberculosis. This example provided by the RATT developers runs smoothly:
Of the reference chromosome Tb_H37Rv    2.66 per cent   has no synteny with the query
Of the query chromosome gi_148821191_ref_NC_009565.1_    3.03 per cent  has no synteny with the reference

...
Please see the file F11.gi_148821191_ref_NC_009565.1_.Report.txt for reporting the errors and the file F11.gi_148821191_ref_NC_009565.1_.Report.gff for reporting the error in Artemis or other sequence viewer.
************************
If you want to start artemis on this replicon:
art F11.gi_148821191_ref_NC_009565.1_.final.embl + F11.gi_148821191_ref_NC_009565.1_.Report.gff  + Query/F11.gi_148821191_ref_NC_009565.1_.Mutations.gff

Then, I tried it on the genome of E. coli and it worked:
$RATT_HOME/start.ratt.sh embl bact1.fa bact1 Species

I used the embl file from Mycobacterium so not much annotation was transferred due to the large genetic differences between Mycobacterium and E. coli.
Then, I tried to transfer the annotation from the reference genome of C. elegans to a new assembly I have generated. The way I generated the embl annotation file is by running seqret on the fasta file of C. elegans and its corresponding annotation file found in GTF format:
seqret -sequence c_elegans.PRJNA13758.WS263.genomic.fa -feature -fformat gff -fopenfile c_elegans.PRJNA13758.WS263.annotations.gff3 -osformat embl -auto -outseq c_elegans.PRJNA13758.WS263.embl

I got the following error:
RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

ERROR: Reference input does not match delta file
Assertion failed: (Len - Coord + 1 > 0), function revC, file show-coords.cc, line 1118.
/Users/user/Documents/project/output/software/ratt-code/start.ratt.sh: line 73: 85325 Abort trap: 6           show-coords -clHT $name.delta > $name.coords
Assertion failed: (Len - Coord + 1 > 0), function revC, file show-coords.cc, line 1118.
/Users/user/Documents/project/output/software/ratt-code/start.ratt.sh: line 73: 85326 Abort trap: 6           show-coords -clHT $name.filter.delta > $name.filter.coords
Use of uninitialized value $stats in print at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 82.
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on c_elegans.PRJNA13758.WS263
Use of uninitialized value in split at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 1140.
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff
work on II
If you want to start artemis on this replicon:
art scaffolds.II.final.embl + scaffolds.II.Report.gff  + Query/scaffolds.II.Mutations.gff
work on III
If you want to start artemis on this replicon:
art scaffolds.III.final.embl + scaffolds.III.Report.gff  + Query/scaffolds.III.Mutations.gff
work on IV
If you want to start artemis on this replicon:
art scaffolds.IV.final.embl + scaffolds.IV.Report.gff  + Query/scaffolds.IV.Mutations.gff
work on V
If you want to start artemis on this replicon:
art scaffolds.V.final.embl + scaffolds.V.Report.gff  + Query/scaffolds.V.Mutations.gff
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

I downloaded the FASTA file here and the GTF file here.
I am working on Mac OSX.
Update 6th August 2018
I have also tried an alternative way of getting a C. elegans annotation in EMBL format. I downloaded it from the Ensembl FTP website:
wget -m ftp://ftp.ensembl.org/pub/release-91/embl/caenorhabditis_elegans/ --directory-prefix embl/
Then, I cat-ed the .dat files (one for each chromosome) into one ensembl.embl file and used that for the RATT command:
$RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

ERROR: Could not parse delta file, nucmer.scaffolds.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.filter.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.delta
error no: 402
ERROR: Could not parse delta file, nucmer.scaffolds.filter.delta
error no: 402
Use of uninitialized value $stats in print at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 82.
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on ensembl
Use of uninitialized value in split at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 1140.
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff

...
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

Indeed, the nucmer delta file is empty. If that is the only problem, maybe I can run nucmer separately but then how can I run RATT without running nucmer?

The developer of RATT told me to not concatenate the EMBL files because there is no need to do so. So, here is the output of RATT on the .embl files, one for each chromosome.
$RATT_HOME/start.ratt.sh embl scaffolds.fa scaffolds Strain

Of the reference chromosome Caenorhabditis_elegans.WBcel235.91.chromosome.I 13.76 per cent  has no synteny with the query

...
Of the query chromosome X    15.29 per cent has no synteny with the reference
Nucmer is done. Now transfer the annotation.
perl /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl embl nucmer.scaffolds.snp nucmer.scaffolds.filter.coords scaffolds
working on Caenorhabditis_elegans.WBcel235.91.chromosome.I
Argument "11123466)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2547.
Argument "11123466)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2555.
Argument "5568506)" isn't numeric in sort at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 985, <F> line 2627.
1474993 1480045 2526 1477519
14358883 14373963 7540 14366423
3204189 3216229 6020 3210209
3809476 3812106 1315 3810791
10232464 10233260 398 10232862
10491944 10506864 7460 10499404
11020347 11023081 1367 11021714
318712 361504 21396 340108
working on Caenorhabditis_elegans.WBcel235.91.chromosome.MtDNA
Can't use an undefined value as an ARRAY reference at /Users/user/Documents/project/output/software/ratt-code/main.ratt.pl line 229.
6 Sequences where generated
Done.
Nucmer is done. Now Correct the annotation for chromosome .
work on I
If you want to start artemis on this replicon:
art scaffolds.I.final.embl + scaffolds.I.Report.gff  + Query/scaffolds.I.Mutations.gff

...
work on X
If you want to start artemis on this replicon:
art scaffolds.X.final.embl + scaffolds.X.Report.gff  + Query/scaffolds.X.Mutations.gff
rm: scaffolds*embl.tmp.BBA.embl: No such file or directory

It is better than before because I get a reasonable delta file from nucmer. But the error Argument "11123466)" isn't numeric is a problem. It seems that is should be easy enough to solve this. It looks to me like the coordinates from the EMBL file are not parsed properly and the code stumbles on the coordiate followed by a ).

The seqret didn't give any error?  It seems like the expectation of ratt and the output of seqret differ. Maybe it needs a specific version or the problem is only on the RATT program

seqret did not give an error.


text:
Can I create a CRAM file with a relative reference path?
I’m trying to create a CRAM file that stores its path to the FASTA reference as a relative path, rather than an absolute path, so that I can move the files around. Unfortunately I can’t get this to work; I was expecting the following to work:
⟩⟩⟩ samtools view -C -T ../reference/ref.fa -o output.cram input.bam

However, the resulting file contains an absolute path in its header:
⟩⟩⟩ samtools view -H output.cram
…
@SQ     SN:1    LN:249250621    M5:hash     UR:/absolute/path/to/data/mapped/../reference/ref.fa
…

As a result, I am unable to open the file via a different path mount that results in different absolute paths, and I can’t move the file (+ its reference) around, or to different machines.
I know that I could set the REF_PATH environment variable or specify -T when reading the file but I would like to avoid this (the result file needs to be readable by IGV, launched by users who don’t know how to set environment variables).
Is there a way of creating a CRAM file that stores a relative path to its reference?

Can you reheader the file? I'm curious if stripping the absolute path works that way.

BTW, ideally the hash should suffice to looking up the reference (assuming IGV supports this).

I have a feeling this should be raised on hts-specs then. What you’re trying to do should really be supported.

This would seem like a rather glaring oversight which won't help with the adoption of CRAM at all, as it effectively makes CRAM unportable which somewhat limits it's utility for an archive format.

The primary way of finding the right reference sequences for a CRAM file is via the M5 hash, which is far more useful in an archival context than a file path on someone else's potentially long-gone filesystem.
If there's an issue to be raised, it's against htslib which unconditionally makes the file path absolute when adding UR from a -T argument. You would indeed be able to set it arbitrarily when reheadering.

In your CRAM file, make all of the reference sequences present in the ./ directory. Then make a bash script to make relative links to your reference for whatever mount you are using at the time - ie ln -s TARGET LINK_NAME.
Now your CRAM is portable and you just have to tailor the local working dir to it.

source: Stack Exchange

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

sort: votes


text:
How do I filter a GFF file by gene type?
I have a file with following information:
chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    141474  173862  .       -       .       gene_id "ENSG00000241860.6"; gene_type "processed_transcript"; gene_name "AL627309.5"; level 2; tag "ncRNA_host"; havana_gene "OTTHUMG00000002480.3";
chr1    ENSEMBL gene    157784  157887  .       -       .       gene_id "ENSG00000222623.1"; gene_type "snRNA"; gene_name "RNU6-1100P"; level 3;
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

In the above file there is a gene_type with different names like processed_pseudogene, lincRNA, snRNA, unprocessed_pseudogene, antisense_RNA, processed_transcript.
I want to create a new file based on gene_type. Among all the gene_types I want the new file to be with certain names of gene_type like only processed_pseudogene, lincRNA, unprocessed_pseudogene, antisense_RNA.
The new file should look like below:
chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

your algorithm to create a file is not clear. please clear your question

Don't you have some bioinformatics tool that filters GFF files on gene type?

No, I dont have any. Are there any tools for filtering on gene_type? Could you please recommend me some.

@Kusalananda They mention that they want to keep "certain names of gene_type like only processed_pseudogene, lincRNA, unprocessed_pseudogene, antisense_RNA"

@Sparhawk Ah, they edited. I missed that. I still think it should be migrated though. I'm fairly certain there are standard bioinformatics tools for these types of operation on GFF files.

@Kusalananda I'm not 100% sure there are, because that 9th field is often just a unregulated dumping ground for "other" characteristics. I'm happy to proven wrong though… it would be most helpful!

@Sparhawk There are mandatory and optional fields and there are controlled vocabularies for that 9th field, so it nat just "anything goes".

grep -E 'gene_type "(processed_pseudogene|lincRNA|unprocessed_pseudogene|antisense_RNA)";' input_file.txt

Explanation

grep -E: use extended regular expression so you don't need to escape things.
This kind of expression (foo|bar|baz) allows grep to include either of these strings. Hence, the total expression will match gene_type "processed_pseudogene";', gene_type "lincRNA";', etc.

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


text:
Installing a R package from git
I am trying to install STEMNET R package but a permanent error stops me, what can I do please?
> > library(devtools)
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=TRUE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=FALSE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > library(git2r)
>
> Attaching package: ‘git2r’
>
> The following objects are masked from ‘package:Biobase’:
>
>     content, notes
>
> The following object is masked from ‘package:igraph’:
>
>     tree
>
> >
> > libgit2_features function ()  {
>     .Call(git2r_libgit2_features) } <bytecode: 0x1d276d918> <environment: namespace:git2r>
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=FALSE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > install_git("c97sr/idd") Downloading git repo c97sr/idd Installation failed: Error in 'git2r_clone': unsupported URL protocol
>
> >
> > install_git() Error in lapply(url, git_remote, subdir = subdir, branch = branch, credentials = credentials) :    argument "url" is
> missing, with no default
> >

I wonder if you have an issue with your SSL library or the server was down, I've just installed the STEMNET package and all dependencies on my mac without problem I'm using R 3.5.1 and High Sierra

UPDATE 2018/04/22: Some users had problem with this command, a way to
  install that works for everyone is to download the tar.gz archive,
  unpack to a directory of your choice (e.g. /path/to) and inside R run
  devtools::install_local("/path/to/stemnet_dir", build_vignettes=T)

https://git.embl.de/velten/STEMNET

Sorry but my problem is this I don't know from where download the  tar.gz archive, do you know?

https://git.embl.de/velten/STEMNET/-/archive/master/STEMNET-master.tar.gz


text:
Are gene names same across species?
I have a bunch of gene names of Apis mellifera (specifically 194). I used these gene names as an input on STRING database to create a network for Drosophila melanogaster. 54 of those genes were also present in D. melanogaster and I received a second network. Is this a credible method of creating homologous gene network. Or should I use the sequences of the 194 genes and BLAST it to the geneset of the D. melanogaster and do it that way?

Wonder why the negative vote. Its a credible question when you see all these papers saying the same genes are present in different species with the same names.

You could try using BioMart in Ensembl Metazoa. Filter by your list of gene names in honeybee, then get the fly homologues as attributes. The homologues are calculated by sequence comparison and clustering.

Gene names are not consistent across species. Most of the time genes are automatically annotated and the algorithms try to give names that are similar to their matches in other organisms, but not always. It's also more likely that genes will be automatically named/annotated based on their sequence, and manually named/annotated based on their function.
The further away you get from Homo sapiens & Mus musculus (and maybe Caenorhabditis elegans, Arabidopsis thaliana & Danio rerio), the less likely it is to find consistently named genes.

Given the edits to the question i.e. honey bee input and expecting idential names in fruit fly homologs, I suspect the optimal approach for the OP would be to BLAST the 194 sequences against fly to work out the appropriate fly gene.

https://omabrowser.org/oma/home/ is useful for finding homologous genes in other organisms, but it looks like it lacks coverage for honey bee.


text:
Hemoglobin subunits genes in scRNA-seq
In one scRNA-seq sample I encountered the genes: Hbb-bs, Hba-a1 and Hba-a2. These genes appear on top of the list of the highest expressed genes but the 75% percentile of cells have very low counts. Therefore the expression is highly skewed.

Reading from gene cards, they appear to be coding hemoglobin subunits and involved in oxygen transport. I think they may need to be removed before normalization, also because the analysis I'm running does not consider oxygen transport.
Has anyone encountered this situation and what do I do?

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

They may be B cells, did you analysis if they have b cell marker?


text:
How I can reproduce this heat map
I have single cell RNA-seq on 8 time points. I have clustered each time point to 2-3 clusters by seurat. I also have a list of differentially expressed genes between all my 1562 cells in 8 time points. I know how to plot this list of genes in heat map on each cluster of cells or a big heat map of my 1562 cells, but I have also sub-clusters for each cluster (time point).
Could you please help me to figure out how can I have a heat map with sub-clusters for each time point on which I am showing a list of differentially expressed genes?
Something like  the picture below where each day has been clustered to 2-3 sub-population then they plotted common differentially expressed genes between day points
Edit
Lets say 500 genes commonly are being differentially expressed in 8 time point. time point 1 has 3 sub-clusters of cells and time points 2 has 2 sub-clusters of cells. time point 3 and 4 both have 2 sub-clusters of cells. I want to heat map these 500 genes in cells in time point 1 (3 sub-clusters), time point 2 (2 sub-clusters), time point 23 (2 sub-clusters) and time point 4 (2 sub-clusters) in one unified heat map on which cells in each subclusteres (a,b and c) are in column and genes in rows but I don't have any idea how to do that. column annotation function in ComplexHeatmap likely works in this regard But I don't know how to define cells inside each  sub cluster in different colours
I can heat map of 500 genes only on time points easily by seurat but without underlying sub clusters
> DoHeatmap(object = SubsetData(object = object, max.cells.per.ident =
> 100), genes.use = 50 genes, slim.col.label = TRUE, group.label.rot =
> TRUE)

Figure 2 in this article.

Some sample expression data to plot in said heatmap will be of use to encourage answers to your question.

What have you tried to do? Do you have each heatmap for each day and with the same order of genes?

I have independently clustered cells in each time point by seurat. I also obtained differentially expressed genes in each time point separately so I have common genes differentially expressed genes between each pairs of neighbour time points.

In case it is still helpful as the post is rather old, the code below would generate a heatmap with annotations thanks to the ComplexHeatmap package.
But before that I would like to stress out the importance of example data, I bet a lot of experts could not help you just because of they did not have single cell data at their disposal.
I am using some scRNA-seq data stored within a Seurat object. The expression data come from the @scale.data slot and the annotations are from the @meta.data slot. I am using a few general cell type markers as genes and 3 columns from the metadata slot (cell type, cell subtype, sample region). When using the ComplexHeatmap package, the trick is to have the same order in the annotation and actual data tables and this is how you should control what (in your case single cells) is plotted in what order.
library(Seurat)
library(dplyr)
library(ComplexHeatmap)
library(circlize) # for the colorRamp2() function

seurat_object <- readRDS(...)

# Annotations

meta_data <- seurat_object@meta.data

# order of annotations/colors are defined here
ordered_meta_data <- meta_data[order(meta_data$cell_type), ]

# OPTIONAL: YOU CAN PICK COLORS FOR EACH LEVEL OF ANNOTATION
# HERE I PROVIDE COLORS FOR ONLY TWO OF THE THREE LEVELS, COLORS FOR THE
# REMAINING LEVEL IS TAKEN CARE OF BY THE PACKAGE
annotation_colors <- list("cell_type"= c("Epithelial" = "red",
                                         "Immune" = "orange",
                                         "Other" = "gray",
                                         "Stroma" = "green"),
                          "sample_region" = c("B" = "blue",
                                              "C" = "yellow",
                                              "N" = "maroon",
                                              "T" = "darkblue"))
ha = HeatmapAnnotation(df = ordered_meta_data,
                       show_annotation_name = TRUE,
                       col = annotation_colors)

# Expression data

genes_to_use <- c("EPCAM",
                  "CD3D",
                  "NCR1",
                  "CD79A",
                  "LYZ",
                  "COL1A1",
                  "CLDN5",
                  "S100B",
                  "KIT")

seurat_object <- ScaleData(seurat_object,
                           genes.use = genes_to_use)

my_data <- seurat_object@scale.data

# COLUMN ORDER OF THE EXPRESSION DATA SHOULD MATCH THE ROW ORDER OF THE
# ANNOTATION TABLE
my_data <- my_data[, rownames(ordered_meta_data)]

# Heatmap

col_fun = colorRamp2(c(-2, 0, 2), c("blue", "white", "red"))
png(filename = "example_heatmap.png",
    width = 1000,
    height = 1000)
Heatmap(
  my_data,
  col = col_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  column_order = NULL,
  show_row_dend = FALSE,
  show_column_dend = FALSE,
  show_row_names = TRUE,
  show_column_names = FALSE,
  use_raster = TRUE,
  raster_device = c("png"),
  bottom_annotation = NULL,
  top_annotation = ha
)
dev.off()


text:
Where to get '--fldMean' and '--fldSD' for single-end Salmon run
I need to run single-end bulk-RNAseq salmon - based alignment and I am not sure what I should supply as --fldMean and --fldSD which seem to be crucial for the single-end experiment alignment. Can they be somehow derived from the FASTQ files that I have, or should I be provided with them by someone who did the experiment?

From looking at the salmon documentation here, it is correct that both of these are essential for single end experiments:

Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads

From the discussion here it would appear that you can't estimate these parameters at all via the FASTQ files alone.
The fragment length refers to the physical molecule of D/RNA captured in the library prep stage which you have sequenced (in part) on your sequencer.  Thus ideally this sort of data should come from a bioanalyzer trace from whoever performed the library prep or similar QC step, this should ideally show the mean fragment length and standard deviation.


text:
Sampling haplotypes
I am trying to simulate different genome of peoples, I have data (VCF files) of various genes from the 1000K Gene project.
I want to simulate different whole genomes i.e generate a new population by combining real haplotypes I have. I am wondering what is the best way to approach the problem. That is what is an efficient method of generating realistic (Not just randomly select parts and recombining them) new genotypes based on the real genotypes I already have.
I am using Bioconductor packages VariantAnnotation to read and manipulate the VCF files and TxDb.Hsapiens.UCSC.hg19.knownGene to determine the Genes positions.
The data looks as the following:
> gene58@rowRanges
GRanges object with 91 ranges and 1 metadata column:
              seqnames    ranges strand | paramRangeID
                 <Rle> <IRanges>  <Rle> |     <factor>
  rs551585351        1 229566998      * |         <NA>
  rs528384854        1 229567027      * |         <NA>
  rs542093083        1 229567063      * |         <NA>
  rs561849701        1 229567128      * |         <NA>
  rs531042647        1 229567160      * |         <NA>
          ...      ...       ...    ... .          ...
  rs565479298        1 229569784      * |         <NA>
  rs572772527        1 229569785      * |         <NA>
     rs605430        1 229569803      * |         <NA>
     rs605428        1 229569804      * |         <NA>
  rs368699658        1 229569810      * |         <NA>
  -------
  seqinfo: 86 sequences from "hg19" genome

The genotype matrix of the gene:
> gene.mat[1:5, 1:5]
            HG00867 HG02371 HG00759 HG00766 HG00844
rs551585351 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs528384854 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs542093083 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs561849701 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs531042647 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"

For simple haplotype simulation of unrelated individuals, there is Montana's venerable HapSim. The model used in that tool is simple and may be good enough for your purposes. You will need to acquire a recombination rate map from for example the HapMap consortium or from 1000G. A cursory google search also revealed this tool and walkthrough that may be suited to your problem: https://adimitromanolakis.github.io/sim1000G/inst/doc/SimulatingFamilyData.html
There are lots of other tools with more advanced models, but it sounds like those would be overkill here.

The data in your matrix appears to be from a VCF file (For reference, see https://github.com/samtools/hts-specs/blob/master/VCFv4.3.pdf).
In the VCF specification, a heterozygous genotype can be specified as 0|1. The zero and one indicate the allele numbers (0=ref, 1=first alt) that make up the genotype. The "pipe" (|) is used to indicate that the alleles are in phase with the previous entry in the file. Therefore, if you observe data like this:
"0|1"    "0|0"
"0|1"    "1|1"
"0|1"    "1|1"

Then this means that the first sample has the haplotypes [0,0,0] and [1,1,1]; and that the second sample has the haplotypes [0,1,1] and [0,1,1].

Thanks, but its not what I meant, I wanted to know how to recombine the haplotypes in order to generate new population out of it. Is there a place where recombination is more likely occur? Is there a program where I can input real haplotypes and it will generate a population?

Oh! Well there certainly are programs for that. Any chance you could rephrase your post to contain that question?

Hopes its better, let me know if I left anything unclear.

Ensembl have haplotypes from 1000 Genomes available by transcript.


text:
Changing the colour of each cell in tSNE plot
I have plotted a tSNE plot of my 1643 cells from 9 time points by seurat like below as 9 clusters. But, you know I should not expected each cluster of cells contains only cells from one distinct time point. For instance, cluster 2 includes cells from time point 16, 14 and even few cells from time point 12. My question is, how I can change the colour of cells here so cells from each time point has a distinct colour does not matter some cells are being clustered somewhere else because this is unavoidable. I have named my cells already so I know cells from each time point but I want to colour cells to nine colours based on the time point they have been sequenced not the cluster they are being grouped by seurat.
I first created a file in one column is name of my cells and in second colour the corresponding time point
Like so
    time
s0.1    h0
s0.2    h0
s0.3    h0
s0.4    h0
s0.5    h0
s0.6    h0
s0.7    h0
s0.8    h0
s0.9    h0
s0.10   h0
s0.11   h0
s0.12   h0
s0.13   h0
s0.14   h0
s0.15   h0
s0.16   h0
s0.17   h0
s0.18   h0
s0.19   h0
s0.20   h0
s0.21   h0
s0.22   h0
s0.23   h0
s0.24   h0
s0.25   h0

I read my file
time <- read.table("time.txt",header = TRUE)
seurat@meta.data$time <- time
TSNEPlot(object = seurat, group.by = "time")

Sorry in my code colour of tSNE changed also the position of cells changed that I don't want that. I want to change the color of cells but they are still in the same cluster
says that
> TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
  'to' must be a finite number

when I changed names(identity)=sampleInfo$cellName to names(identity)=src$cellName
that worked but again the position of cells changed that I don't want that :(

If you can extract the underlying data (as this is essentially an xy scatter plot) you can then share it here.  Your question then essentially becomes a base R / ggplot2 problem which should enable others to help you better than the screen shot alone.

From an experimental interpretation point of view if your cells are failing to cluster by time point then clearly there is something else going on with your data, again more information would be useful, additionally if you want help with plotting we'd need to know which time point each data point is from.  Thus some exemplar data is required for us to be able to help.

Please  no negative point :( If I lose my reputation I can not ask question anymore

I found the solution   if this is my time table data cc <- AddMetaData(object = n, metadata = time)
cc <- SetAllIdent(object = cc, id = "time")
TSNEPlot(object = cc)
This only colours the cells based on the time point they are

Seurat puts the label in the tSNE plot according to the @ident slot of the Seurat object.
By default, if you do the tSNE without computing the clusters and you have the correct metadata in the object, the labels should be pointing to your timepoints not to the clusters.
If it is not the case, what you have to do its to change this slot to reflect to which timepoint belong each cell you are plotting.
For instance if your Seurat object is called src:
identity=as.factor(sampleInfo$timepoints)
names(identity)=sampleInfo$cellName
src@ident=identity

And then plot the tSNE. Something like this should work after adapting it to your metadata.

Sorry what is sampleInfo?

What you call time in your example. cellName could be your s0.1, etc and timepoints your h0, etc

Sorry says that  > TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
'to' must be a finite number

Sorry I can nit understand. I have the name of my cells and time points as I showed in my question, so sampleInfo is name of my time contains cell names and time points?

Yes, it is supposed to contain similar information. I put this name randomly. In your case it would be your time data frame

Sorry but says that    > TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
'to' must be a finite number


text:
Reproducing GTEx transcriptome analysis
I am willing to reproduce part of the analysis from "The human transcriptome across tissues and individuals" (Melé et all, 2015).
I downloaded GTEx v6 FPKM data in txt format from GTEx portal. I want to mimic the MDS plot from Supplementary Figure 3 (or get the equivalent PCA).
In my process I did:

Load the FPKC matrix
Remove those transcrits not having a level of 0.1 in minim 10 samples
log10 transform after adding 1 to the matrix (to avoid 0s transformation)
Quantile normalize the matrix
Compute the distance matrix between samples
MDS and plot

But my plot has nothing to to with the one in the paper.

Can anyone see any error in my process o see what I am doing that is totally different from GTEx process pipeline?

Extra info
library( preprocessCore )
library( smacof )

# Load gene expression matrix from file
gen <- read.delim( "data/v6/GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct", skip = 2, heade = TRUE )

## First two columns corresponds to "Transcript name" and "Gene name"
gen_ann <- gen[ , c(1,2) ]
gen_clean <- as.matrix( gen[ , -seq( 2 ) ] )

# Load map file
sample_map <- read.delim( "data/v6/GTEx_Data_V6_Annotations_SampleAttributesDS.txt" )
sample_map$SAMPID <- gsub( "-", "\\.", sample_map$SAMPID )
rownames( sample_map ) <- sample_map$SAMPID

# Match data-sets
common_samples <- intersect(
    colnames( gen_clean ), rownames( sample_map )
)
gen_clean <- gen_clean[ , common_samples ]
sample_map <- sample_map[ common_samples, ]

# Select tissues
tissue_to_get <- c( "muscle", "brain", "heart", "blood", "pituitary", "adipose tissue", "nerve", "lung", "breast", "blood vessel", "thyroid", "liver", "skin", "testis" )

samples_to_use <- sample_map$SAMPID[
     tolower( as.character( sample_map$SMTS ) ) %in% tissue_to_get
]

# Filter gene expression matrix according to GTEx
trs_to_use <- apply( X = gen_clean, MARGIN = 1, FUN = function( row, exp_th, min_sam ) {
  sum( row > exp_th ) >= min_sam
}, exp_th = 0.1, min_sam = 10 )
gen_rd <- gen_clean[ trs_to_use, samples_to_use  ]

# log10 transformation and normalization
gen_norm <- normalize.quantiles( log10( gen_rd[ , list_samples ] + 1 ) )

# MDS
d <- dist( t( gen_norm ) )
mds <- mds( d, type = "ratio", ndim = 3 )
[ ... creation of color map ... ]
plot( x = mds$conf[ , 1 ], y = mds$conf[ , 3 ], col = clr2, pch = 20 )

Welcome to the site Carles! Are you comparing your figure with the figure 3 A)  or Figure 1 B) (left)? The step 4 you mention is not used in the article, also they filter without any consideration of the number of samples it is present

Hi @Llopis, thanks for the welcome! I am comparing my MDS with the one from Fig S3 (Multi dimensional scaling (MDS) of GTEx samples based on gene expression). I agree that step 4 (quantile normalization) is not described but I assume they apply some kind of normalization. And they apply quantile normalization in the array set, that is the reason I applied it.

Ok, then my other suggestion is to do the distance metric as they do, 1 - Pearson correlation as they did, maybe you get closer results.


text:
Interpretation of the violin plots from sc-RNA-seq
I'm confused about the meaning of the black dots and the red shape in the violin plots from the seurat tutorial:

It would help if the reference, or legend to this figure was included in the question.

It's from the Seurat tutorial  https://satijalab.org/seurat/pbmc3k_tutorial.html

The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.
The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.

is it normal that you can only see the dot but not the red shape after you doing the Vlnplot? My data shows that problem after I doing the gene in cluster, so I was confuse whether it is a problem or not?

If you see just a dot, it probably means you have one outlier. If you look closely, you will probably notice the rest of the dots at 0 (so they look like a line).

The plot includes the data points that were used to generate it, with jitter on the x axis so that you can see them better. In red you see the actual violin plot, a vertical (symmetrical) plot of the distribution/density of the black data points.


text:
Find a map/correspondence between two versions of a genome
I am working with two versions of the C. elegans genome. I am finding interesting regions (specifically, tRNA genes) in version 1 and then I would like to know if version 2 also has a tRNA gene in that region. For this, I need to know the correspondence between coordinates in genome 1 and coordinates in genome 2. I think a whole-genome alignment programme might help me, e. g. nucmer or mauve. I have run nucmer on the two genomes and got a delta file but I don't know how to interpret it. I have run show-coords on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:
/Users/user/Documents/project/output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS245.genomic.fa /Users/user/Documents/project/output/genome/celegans/strain/seq/bristol/scaffolds.fa
NUCMER

    [S1]     [E1]  |     [S2]     [E2]  |  [LEN 1]  [LEN 2]  |  [% IDY]  |  [LEN R]  [LEN Q]  |  [COV R]  [COV Q]  | [TAGS]
===============================================================================================================================
       1      432  |      498      929  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I
       1      432  |     1002     1433  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I

Does this mean that chromosome I starting at 1 and ending at 432 of version 1 corresponds to chromosome I starting at 498 and ending at 929 and that there is 100% identity between these two 432-long sequences? In this case, it's easy to find a correspondence between bases in version 1 and bases in version 2. However, some lines do not align sequences of the same length:
   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I

Also, in this sequence position 1 of chromosome I is aligned again? So is it aligning to more than one position in version 2 of the genome? Also one sequence is 237139 bases long and the other 237160 bases long. How do I know the internal structure?
Can anybody guide me, please?
I've also aligned the two genomes using Mauve:
progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve > scaffolds.txt

I've tried to import the .mauve file into Python using Biopython but I got an error:
alignment = SeqIO.read('scaffolds.mauve', "mauve")

Traceback (most recent call last):
  File "code/mauve.py", line 4, in <module>
    alignment = SeqIO.read('scaffolds.mauve', "mauve")
  File "/Users/user/anaconda/lib/python3.6/site-packages/Bio/SeqIO/__init__.py", line 712, in read
    raise ValueError("More than one record found in handle")
ValueError: More than one record found in handle

Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.
I am working on Mac OSX.
I would appreciate your help in understanding the delta file, resolving the Python error or any other support in finding a function between coordinates in version 1 of the genome and coordinates in version 2 of the genome.

SeqIO.read is for a single record, try SeqIO.parse

NCBI has  REMAP page/tool where you can feed it a file with coordinates of sequence features from one build version, and ask for those features mapped to a different build.  I have only used it to go forward, not backword.  Also, each version of the C. elegans genome used to have a liftover directory at the Sanger ftp site that showed the version changes, as I recall. Check out the WormBase site.

I don't see where I can feed my own file to the Remap service of NCBI: https://www.ncbi.nlm.nih.gov/genome/tools/remap. Reply from NCBI to my request: "The remap service relies on assembly to assembly alignment data to do this remaping/conversion. These datasets came from NCBI's eukaryotic genome assembly workflow. For outside genomic assembly that are not deposited to NCBI, we have no access, so we will not be able to provide mapping between that and those available from the NCBI collection."

The file that you need is a bed file, or a GFF3 file, or similar, that contains the start and end coordinates for all of your tRNA gene regions from genome Version 1.  You select and upload that file and tell REMAP the version/build using one drop down menu. Then in another drop down menu you select the genome version/build for Version 2. Click on the button and in a little while you will have some files to download showing the coordinates for those gene regions in Version 1, and how they map to Version 2.

@mdperry the versions I am using are personal. I have reassembled the genomes using my own data so I cannot choose a version available on NCBI because I have not deposited my versions yet.

Hmmmm, okay; that was not clear to me from your original post.  I see that one of the files you mention has ‘WS245’ in the name and I assumed that was a specific official build of the C. elegans genome, and I erroneously concluded that the other ‘version’ (you mention version 1 and version 2) was a different specific official genome build, and my comments were for that particular situation.  If you are using your own, homemade versions wont you eventually have to map the genes of interest back to one of the official genome builds (i.e., prior to publication)?

That's a good question. For the time being, I use my own build because I think it's better than WormBase's. But soon, they'll publish a better version and I'll use that one for the publication. However, I'll still have a situation of one official version + personal version. As far as I understand, REMAP only works with official versions.

A relatively convenient method would be to create a .chain file, as used by tools like liftOver and crossmap (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.

I am following all the steps. I am stuck here: faToTwoBit ../../genomes/genbank/GCA_000004515.3_Glycine_max_v2.0_genomic.fna.gz GCA_000004515.3_Glycine_max_v2.0.2bit
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
Referenced from: /usr/local/bin/faToTwoBit
Reason: image not found

I suggest you install that with conda instead, then you'll get the mysqlclient dependency (presumably). BTW, I expect you need to gunzip the fasta file.

I have the following error: "ssh: connect to host localhost port 22: Connection refused" for the command "time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=pwd
-ooc=pwd/${target}.ooc -fileServer=localhost -localTmp="/dev/shm"
-bigClusterHub=localhost -dbHost=localhost -workhorse=localhost
-target2Bit=pwd/${target}.2bit -targetSizes=pwd/${target}.chrom.sizes
-query2Bit=../${query}/${query}.2bit
-querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1". I don't understand why it's trying to ssh.

I don’t know where ssh is coming from there either.

This post helped with the ssh error: https://stackoverflow.com/questions/17335728/connect-to-host-localhost-port-22-connection-refused#22255174

Ah, it uses ssh to submit commands to itself as an ersatz cluster.


text:
Genebank Submission
I am not sure this kind of question is allowed here. I am asking because some of you might have already experienced this.
I submitted some virus genomes in the Genebank. The genomes are on hold and are not public yet.My PI thinks thank we should rename the isolates. I have sent two emails to the admin asking that I want to rename the isolates. I sent an Excel file containing accession number, current isolate name and new name that I want to assign. It has been almost a month and I haven't heard back. So, I don't know if this is the right approach.
What should I do now ? Should I just resubmit the genomes ? Or is there any other way to do that ?
Thanks

I usually hear back from the GenBank curators within a single day.
That said, I’d still suggest you wait: it’s unlikely that they missed both of your emails (but maybe double-check the address and whether they’re likely to be caught by a spam filter): they’re probably just swamped with other submissions.
Furthermore, it’s holiday season so they might operate at drastically decreased capacity.


text:
How to generate Genbank format out of paired end fastq data?
I've unmapped cleaned paired end sequence data in fastq format of a bacterial genome. I want to get a sequence data in Genbank format in the end. What are the exact steps that I have to follow in galaxy webserver and if in some point necessary, also using linux based apps? I would appreciate a lot of details.

The steps needed to assemble a genome in Galaxy are layed out in the Assembly training material from the Galaxy Training network. The software for this should be available on the main galaxy servers (e.g., usegalaxy.eu). Since you have a prokaryote, I suggest you follow the "Unicycler Assembly" workflow, which is also available on usegalaxy.eu.

I think the reason for your confusion is vocabulary.
The keywords you search for are genome assembly (see @Devon Ryan's answer) a process that is trying to stick together reads to recreate the sequenced genomic sequence. The assembled sequence needs to be annotated if you are interested where are the genes or what they could possibly do. If you will search for tutorials using these keywords you will find plenty.
I am not sure why you need it, but Genbank format contains the genome (the product of genome assembly), metadata (who is the submitter, what species it is, etc) and genome annotation (what is in the region of genome).


text:
RCurl issue when installing SCnorm
I am working on the Rivanna cluster that is of SLURM-type. Because of that I do not have write permissions and need to install the newest R version myself locally, and use it. For it to be successfully built from source, I needed to also build from source libcurl4-openssl-dev, so I downloaded tar.gz file from here:
https://packages.debian.org/ru/jessie/libcurl4-openssl-dev
and did:
mkdir build
./configure --prefix=/scratch/nv4e/curl/build
make
make install

Then, in the same way I built from source R 3.5.1, the latest current version. It was compiled without any errors, but with some warnings. Then, I submitted the R script that uses SCnorm package:
https://bioconductor.org/packages/release/bioc/html/SCnorm.html
Inside of the script I am having:
source("https://bioconductor.org/biocLite.R")
biocLite("SCnorm")
biocLite("SingleCellExperiment")
library(SCnorm)
library(SingleCellExperiment)

Of course, I am submitting it through .slurm file where I have the main command:
Rscript --vanilla scnorm.R

It fails, apparently because SCnorm requires RCurl, but it can not get one, even though in the log file it is trying to access it correctly:

The libcurl library was built from source in the /scratch/nv4e/curl/build folder. The only suspicion that I have here is that maybe its version is wrong and SCnorm needs another libcurl version? Any suggestions would be greatly appreciated.

Update

I tried reinstalling libcurl library from the official website, using the latest version:
https://curl.haxx.se/
But it did not fix the issue. I have also tried using libcurl from Anaconda but it gave the same error. So, neither 7.61.0 (latest), nor 7.29.0 (Anaconda), nor 7.38.0 (Debian link) are working.

This is a bug affecting the curl compile, it's been reported here and also here, I suggest trying to compile from the newer release of the libcurl source rather than the .deb you linked.

I tried reinstalling libcurl from here: https://curl.haxx.se/, and it did not help. I see that the R is using the newly reinstalled version, and giving the same error.

Are you sure it's not hitting the old system version at compile time?

Yes, sure. Cluster support fixed it. The issue was related to gcc version, but I am not completely sure what they did.

Glad it got fixed, although if you find out exactly it would be good to post and answer since this will help others with the same issue.


text:
Parsing gene expression value for a list of genes from GEO
I am looking for the gene expression data in GEO GSE15543 with the following keywords,
(Pancreatic Beta Cells) AND "Homo sapiens"[porgn:__txid9606]

I would like to parse the gene expression value of a list of genes from all the samples reported in any particular study.I understand there are two different formats for downloading datasets, SOFT and MINiML.
I have downloaded the soft.gz file and I am using BioConductor to read the contents from soft file using the instructions given here . The file has the content that I 'm looking for , I have checked the contents on the geo2r webpage. I'm still facing challenge in finding out how to fetch out the gene expression values for a list of genes. Any suggestions?
From the files that I download, I wish to parse the gene expression value and the Sample Title and Sample Description so that I can classify the values based on disease state and control state. Any suggestions on which format can be downloaded for the above-mentioned task?

GEO provides a GUI for this (just click on the "Analyze with GEO2 button"), you can even check/download the R/Bioconductor code used for the analyses.
Basically you will create groups of samples to compare/contrast and R/Bioconductor will run behind the scenes. Here is a nice tutorial.


text:
Path to career in bioinformatics?
I've been completing various courses on bioinformatics through edX and Coursera, and I've made a strong footing in the material.
I'm not sure what to pursue next in bioinformatics that would bring me closer to research in the field or employment.
I don't have bioinformatics at my university and thus am unsure of how to progress.
I have a background in mathematics, chemistry, physics and computer science. So, what would you recommend I do next?

Welcome to the site Rock910. The site is about questions while doing bioinformatics not about the career. I think it would be better suited for [meta]. However I would say that what you need now is to find a job or start doing some analysis in bioinformatics (Although some more background in biology wouldn't hurt). But that depends on what do you like and what do you want to do

Thanks @Llopis, I'll ask it in the meta, I appreciate your comment and advice.

@Rock910 try also https://www.reddit.com/r/bioinformatics/ which as more information about that subject - plenty of old threads with similar questions to yours.

Do you ask on bioinfo career or good bioinfo materials? One belongs to academia.SE and the second here.

The Bioinformatics Stronghold at Rosalind is an excellent resource for building your understanding of many concepts and algorithms you'll see frequently in bioinformatics and genomics.
I'd also say that, if you can make it work, an internship where you can "get your hands dirty" with real problems is a very valuable experience. This option may not be available to everyone due to a variety of factors, but when it is an option it provides a tremendous opportunity to build experience with little risk to the intern or the employer.


text:
Filtering the expression values of a gene from GEO
I'm filtering the gene expression value of a particular gene using Bioconductor in R. Following the instruction given here, I could filter the gene expression value for a gene from the dataset GDS3782.
The output is as follows for the gene of my interest,
GSM524151  "23.82"
GSM524152  "27.49"
GSM524153  "24.26"
GSM524154  "17.98"
GSM524155  "25.87"
GSM524156  "26.16"
GSM524157  "26.02"
GSM524158  "25"
GSM524159  "29.85"
GSM524160  "18.25"
GSM524161  "28.86"
GSM524162  "27.33"
GSM524163  "27.42"
GSM524164  "28.23"
GSM524165  "22.09"
GSM524166  "24.85"
GSM524167  "27.11"
GSM524168  "32.01"
GSM524169  "28.97"
GSM524170  "30.98"

The above data corresponds to disease and control which I get using
Meta(gds3782)$description

I wish to see the description of each value. For example, from the webpage, I have seen that the first ten values correspond to "disease-state" and the last ten correspond to "control". How do I view this classification while processing the data using Bioconductor in R?

Welcome to the site. First, is your matrix a matrix of numbers or it has text too? Second in the Meta(gds3782) is there any column describing this or not? (It should but I don't know how it is named in this case). You need only to locate the column

You're looking for the Columns() function:
> Columns(gds3782)
     sample        disease.state           description
1  GSM524151 non-diabetic control  Value for GSM524151: beta-cells_non-diabetic condition_donor1; src: human beta-cells_non-diabetic
2  GSM524152 non-diabetic control  Value for GSM524152: beta-cells_non-diabetic condition_donor2; src: human beta-cells_non-diabetic

Note that it's probably easier to work with a more common structure like and expressionSet:
eset = GDS2eSet(gds3782)

You'd then use the pData() function to return the same thing.


text:
RNA-Seq type and and optimal fusion detection
There are several popular types of RNA-Seq library prep which are frequently used: total RNA (with and without ribosomal depletion), mRNA-Seq/poly-A, and targeted mRNA-Seq/RNA exome. What would the impact of choice of these various types of RNA-Seq have on downstream fusion detection by programs like: FusionCatcher, STAR-Fusion, Arriba, and Pizzly?
For instance whist these methods have various biases and known issues, in particular rRNA depletion badly effects the accuracy of comparative exon level transcriptomics.  That and RNA degradation is known to affect the efficacy of fusion detection. To date has anyone studied directly the choice of methodology on the performance of fusion detection directly, or can we infer that one method should be better than another for this usage case?
Obviously we can assume that targeted approached will only find partners of known genes, but will the different physical capture techniques introduce biases or shortfalls?


text:
SC-RNA seq percent.mito
Is that possible that the percent.mito level is zero? I have done vlnplot and it shows that its percent.mito level is zero. Does that mean my  mitochondrial gene didn't detect?

It is violin plot. And the command is:

Sorry hua, but I find it difficult to understand. Could you clarify what is "percen.mito" ? It seems like it should be a function or a property of your data but I'm not sure. Also, could you explain what is "VlnPlot" ? It seems like a violin plot, but I'm not sure, perhaps the image didn't get attached (upload just a link and we will add it).

Your genome likely doesn't title mitochondrial genes MT-something, so you'll need to change how mito.genes is created.

Is "percent.mito level" the percent of reads mapping to genes annotated as mitochondrial? Are the mitochondrial genes annotated? Note you have [edit] button you can use to improve your question.

You should put "Seurat" in your question, and add a Seurat tag.

I've just come across a similar issue. It seems the mouse and human mitochondrial genes use mixed capitalizations. See the following PR for details:
https://github.com/Hoohm/dropSeqPipe/pull/43
Bottom line is, if you ignore the case in the grep pattern it might work (assuming you use human/mouse):
mito.genes <- grep("^mt-", rownames(shox2crepos_mtmg@data), value=TRUE, ignore.case=TRUE)

Thank you so much, I have solve this problem after your instruction.

The simplest answer is that the bit of code you copied from the Seurat vignette is not working on the gene names you have.  Look at the code.  If none of your gene names begin with "MT-", none will be identified as mitochondrial.


text:
DiffBind: Model must include count data for contrasts.
I am trying to use DiffBind for my ATACseq experiment going through the vignette for DiffBind. I was able to import the data and make a heatmap with plot(data). But somehow calling data <- dba.contrast(data, categories=DBA_TISSUE) function is not working for me giving a warning message:

Warning message:
  Model must include count data for contrasts.

I tried running beforehand: data <- dba.count(data) but it is giving an error:

Error: Error processing one or more read files. Check warnings().
  In addition: Warning messages:
  1:
  2:
  3:
  4:
  5:
  6:
  7:

My sample_sheet.csvlooks like that:

I just need to compare Dura vs Microglia which is Tissue column. How could this error be fixed? Any suggestions would be greatly appreciated.

Note that you have a BED file in the bamReads column rather than the BAM file that you fed to MACS2...

I know but it is written in the documentation that both bed and bam files are fine. I tried supplying bam files too since I have them, but it does not work either.

Sure it can accept BED files, but they need to be from something like bedtools bamtobed. What was the error when you used BAM files?

Wrong suffix error and just empty warnings as shown at the top. My data quality is poor (low mapping rates), so can it just be data quality issue?

Could be. You might try CSAW instead, it's what we use for differential accessibility with ATAC-seq data.

Can you provide some reproducible data to test and see how to avoid the error?

With the data that I have now it is definitely not working, and the issue is not yet solved either because of the poor quality data, or because something else is going wrong. I will post the answer if I figure out why that was happening. Besides, I am not sure how I can attach 12 files here, or even just 2 for 1 sample.

Ah, I have definitely run into this issue before - it is resolveable, but I cannot for the life of me remember exactly what the issue was. I think it was related to the file paths, but your data quality is likely not the issue.

The problem was that I was using wrong files. I ran a pipeline from here and it generates several folders with many files. The .bam files that should be used are located in the aligned_mm10 folder and are with the extension: .pe.q10.sort.rmdup.bam.

This link leads to a page not found. What was the pipeline that you were using?

I see, it is just a private page, so that is why you might not be able to access it. It is this one: https://github.com/databio look for pepatac. Then, my answer is useless since no one can access it anyway...


text:
Glocal\semi-local\Hybride Globale-Local alignment with Python
I was looking for a simple way to do a glocal alignment.
The case I have is I have a small sequence which should be find in a bigger one, thus typically a glocal alignment.
Also I can not Install anything on the target machine (Python 3.6 with Biopython is already installed), all I can do is pip install.
My question is : are there any glocal alignment implementation out there ?
I was unable to found one and I don't feel proficient enough with Python yet to implement it myself.

@winni2k there is also a newer pairwise aligner in Bio.Align http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc87

@ckin, is there any particular reason why you need to perform a glocal alignment?

@winni2k like I said I have a typically case where you want to do a glocal alignement: I have to align a small sequence against a bigger sequence, and my small should be found within the bigger one.

As suggested in the comments, the built-in BioPython pairwise aligner should work for this. You can find a short sequence match in a longer sequence using a global alignment with the right score parameters (query_end_... = 0). The other scoring parameters are values that I've found to be good for aligning similar sequences.
Example code:
from Bio import Align

aligner = Align.PairwiseAligner()
aligner = Align.PairwiseAligner()
aligner.match_score = 5
aligner.mismatch_score = -9
aligner.mode = 'global'
aligner.target_internal_open_gap_score = -12
aligner.target_internal_extend_gap_score = -3
aligner.target_end_open_gap_score = -12
aligner.target_end_extend_gap_score = -3
aligner.query_internal_open_gap_score = -12
aligner.query_internal_extend_gap_score = -3

# Set the query_end penalties to zero to effectively implement a 'glocal' alignment as described in the question
# These are the default values, so you could skip the next two lines,
#    but this is relevant to the question, so I'm including them in the example.
aligner.query_end_open_gap_score = 0
aligner.query_end_extend_gap_score = 0

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATCCG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACG")

for a in alignments:
    print(a)
print()
print()

alignments = aligner.align("GAACTCGCATCGATGATCCGTGCGCTGAGATCCGCTGACATGATC", "GATACCG")

for a in alignments:
    print(a)


text:
How to plot circular chromosomes with Gviz?
I am using Gviz to plot several tracks. I create a GenomeAxisTrack() and then several AnnotationTrack()s that I feed to plotTracks(). This produces a stacked output of tracks with annotated arrows as I want but the displays are linear. In some cases I have the complete genomes and I know they are circular (like some bacterial and viral chromosomes), how can I display the plots as circular?
I tried feeding to AnnotationTrack(); a GRanges() object with a Seqinfo() object specifying the isCircular=TRUE flag, but this didn't change the resulting plot.
Here is minimal code to create a linear plot, I would like one where the first and last arrow join up in a circle:
library(Gviz)
library(stringr)

track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('+','+','-'), stacking="dense", fill=c('red','blue', 'black'))
png(file=paste("test","png", sep="."))
plotTracks(track)
dev.off()

Note, I am not interested in alternative software recommendations, I would like to solve this with gviz if it is possible. Also asked this here.

I would be amazed if there were a way to do that.

@DevonRyan why? i have use cases, but i'm not very familiar with gviz and its implementation

For track-based plots you would pretty much have to set it up for circular plotting from the very beginning in order for it to work. I'm aware that there are use-cases, but given where gviz was developed and the use-cases there it would be surprising if circular genomes were a consideration during development.

The primary author of Gviz confirmed to me via email that this is not possible:

Gviz will only give you the traditional linear tracks. There are a
  bunch of solutions available for visualizing circular genomes, most of
  which are more interactive in nature and not implemented in the rather
  static R graphics context. GGplot2 has some support for circular axes,
  too. If you want to stay within the R environment I’d start there.
  Ggbio is another option providing a somewhat higher-level interface.


text:
Comparing the gene expression data
I'm trying to compare the gene expression data reported in the studies documented in GEO for a specific gene expressed in a tissue of my interest in Homo sapiens.
I compared the values reported in 4 studies, in two studies the expression value of gene A lies between 15-30, in the other two studies the range lies between 1-4.
I would like to ask for advice on whether it is meaningful to compare the data from different studies. What are the factors that one should keep in mind while making such comparisons?
Any help would be highly appreciated

From the suggestions below, as a first check, I have looked for the platform used for the microarray analysis. The platform varies in each study. I'm trying to compare the probe densities. Unfortunately, I am not able to find this information.Any suggestions?

You won't be able to make comparison in this case.

Could you please explain why?

The two answers here.

There's going to be a batch effect between studies, so unless you have some sort of constant group that you can use for between-study normalization you're going to completely muck up your results from using data from multiple studies.
More specifically, this is what packages like sva and the combat() command in it are used to correct. Note that you should go back to the raw .CEL or .fastq files if you want to process multiple studies, in order to at least ensure that everything is processed in the same way.

Please check my recent comment in response to your answer

It doesn't make sense to compare gene expression data from different studies unless you have a negative spiked-in control such as sequins.
The most obvious thing you will need to do is normlizate the data for various bias such as batch effects, sequencing platforms, sequencing depth, bioinformatic workflows and more ... Your task will be quite complicated as you must justify how your normalization would work for data from different studies.
You should check how the raw data was generated (methods section), convince yourself, and then apply gene quantification yourself from the raw files.

Please check my recent comment in response to your answer


text:
fuse fastq files with multiple records
I'd like to fuse ~50 pairs of corresponding fastq files on per sequence bases. I'd be happy for a solution as a Linux script, R script (if feasible with huge files), or any free software that runs preferably under Debian Linux (otherwise Windows). Aim is preprocessing for a Linux pipeline for NGS data. Ideally both the sequence and quality data are fused.
The aim of this request is the following: there is an existing pipeline for downstream mapping of immunoreceptor sequences (migec/mixcr/vdjtools etc.). Previously we got Illumina reads (read 1, read 2 = reverse), with the index and UMI sequences attached to one of the reads => can be alligned and processed to habe 1 biological useful sequence with UMI and index at one end. NOW we get read 1, read 2 and seperately "read 3" which is just the UMI+Index reads, and which was attached by the core facility to read 1 using software we do not have access to. Since the previous pipeline works well, we prefere to use it in the future, but are open for suggestions how to use the index reads with migec/mixcr/vdjtools otherwise.
Each of the 2 files contains the same number (millions) of short (~300 bp) and corresponding DNA reads:
File_A.fastq:

@record1 / aaa... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

File_B.fastq:

@record1 / ttt... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

The target file should look like:

@record1 / aaa...ttt... / +... / ...
@record2 / ...
@record3 / ...
... (millions)

(of course the aaa...ttt... is just an example)
I was looking for a solution on Stackexchange and Google, but the solutions, I found do not work for me because

manual processing is not an option given the number of sequences
I need to fuse the files on a per sequence basis, i.e. a simple Linux command like cat file1 file2 is not what I need (as far as I understand, this would give me (record1_A record2_A.... record1B record2_B).
R, package shortread, readFastq + writeFastq would be perfect, but cannot handle such large numbers of input lines

Maybe I did not use the right search terms, but in any case I'd be happy for a solution

If you want to join the forward and the reverse sequence, why do you want one sequence to be after the other one? Simply it isn't that way in the cell so most programs don't aim to do that. If these are not forward and reverse sequence, did you try perl, python or awk? Are you open to a solution in those languages?

Thanks - I added an explanation as second paragraph. These are of course not fw/rev reads. I have unfortunately not used python, perl or awk so far (but am open if it comes with an easy explanation)

From my understanding  of the problem (which might well be wrong): If the sequences are in the same order (i.e. lines correspond to each other) one very simple approach would be just pasting the files together with paste -d "" File_A.fastq File_B.fastq > target.fastq. This would of course require some clean up first.

@SebastianMüller: Yes, seq are corresponding and in the same order. And they have the same correspnding names in both files (sorry - I made this more precise in the question), i.e. only the sequence and quality data (i.e .lines 2 + 2n), but not the names (lines 1 + 2n) should be fused. Is there a way to do that, or afterward replace fused names with the original ones by another command? (I am new to Linux, was looking for sth. like this but could not find). If sth. like this works with 10+ mio seq per file, this would be the desired answer!

Based on the comments of the question I thought to put together this answer using GNU coreutils which should exist on pretty much all systems.
The idea is to paste the 2 files together with paste (assuming they are in the same order and each line correspond to the same line in the other file).
Since only the sequence and quality of both files should be kept, the lines starting with @ and + need to be made empty for one of the files. This can be achived with sed, e.g. sed 's/^@.*//' empties all lines starting with @ but keeps the empty line.
To modify the second file on the fly, one can using process substitution <().
Everything strung together should be something like this:
file1=File_A.fastq
file2=File_B.fastq
paste -d "" $file1 <(cat $file2 | sed 's/^@.*//' | sed 's/^+.*//') > target.fastq

I haven't tested it thoroughly, but could be a good starting point.
If you have a lot of files to paste, you might however want to process them first (not on the fly) using:
for file in $(ls *.fastq); do
  cat $file | sed 's/^@.*//' | sed 's/^+.*//' > ${file%.fastq}_cleaned.fastq
done

This creates File_A_mod.fastq from File_A.fastq etc.
Best to move them in a separate directory (make sure the first one is the original instead of the _cleaned one since it serves as the base).
You should have the following in the folder: File_A.fastq File_B_cleaned.fastq File_C_cleaned.fastq ....
Those can be than pasted with paste -d "" *.fastq > target.fastq.

What you actually want to do is overlap the reads, for which you can use tools like FLASH or BBMerge (I've had better luck with FLASH).

Thanks, the 2 sequences I want to process are NOT fw/rev or any other overlapping sequences (see new explanation above)

Ah, but at least mixcr will accept the two fastq files as they already are.

Yes, but somehow I need to keep the UMI/index from the 3rd file... If there is a way to do this with  mixcr&co, that would be perfect, but I am not aware of that possibility.

The samples should have been split according to index (I assume this is an Illumina sample index) and the UMI is typically added to the read name (I or someone else can produce a few lines of code to do that for you).

Thanks a lot. Yes, it is Ilumina. Unfortunately sequences are not split for index, and index +UMI are not included in filenames (will add examples as soon as I get faster internet), but Migec can list each index in a barcode file and demultiplex according to index, if supplied with 2 (not 3) files, with UMI+index added to e.g. read 1... (used to work perfectly with this pipeliine as long as we had index+UMI added to read 1)...

Make your life easier and ask your sequencing provider to just give you the files in a standard format (e.g., split by index and with the UMIs in the header). Darn near every tool expects exactly that.


text:
RNA-Seq: clustering/treatment of genes with low expression
I have some RNA-Seq data from leukaemia patients. I want to do unsupervised clustering on them with some other published leukaemia RNA-Seq data and see how they cluster. There are a few problems I encountered while doing this.

I read mix messages of whether using log2(TPM) or rlog(counts) (e.g. Deseq2 rlog or limma voom transformation) for clustering. Which one should I choose?
If I am filtering out genes with low counts, should I do it prior normalisation to library size or after?
I tried using filterByExpr from edgeR for filtering but it removed many of the genes that are only expressed by subtypes with small sample number. If I am to filter the counts "manually", is there a recommendation of how the threshold should be set?
If I should use log2(TPM) for unsupervised clustering, can I treat low TPM (e.g. <2) as 0? As this paper suggested TPM < 2 basically means no expression, is it alright to treat them as 0?

rlog(normalized counts) is going to be more robust than log2(TPM), so use it instead.
Do it afterward, keeping the low counts can be helpful for library-size normalization.
This is very much more of an art than a science, though have a look at the kOverA() function from the genefilter package. That will give you a bit more fine-grained control on filtering.
You'd be better off using something like zFPKMs if you want to filter by "expressed or not", since there's no really great definition of that.

Thanks Devon! Your answers clear things up for me.  May I also ask at what point is transformed TPM/FPKM/other metrics of same sort be useful? I think this is also something bothering some wet lab biologists who just start doing bioinformatics.

FPKM/RPKM should essentially never be used. TPMs would work fine here as well, but more so if you make them from normalized counts. The general thing with these transformed metrics is that they try to take gene/transcript length into account. That's fine for some things but it'll change how you might interpret clustering results (at least if you don't scale things).


text:
Minimap2 -ax map-pb doesn't output tlen field
I have used minimap2 to map some pacbio reads to a reference genome. I would like to know the "insert size" (true length of sequence) relative to the reference. More specifically, I want to know the first base position in the reference where the read maps and the last base position in the reference where the mapping ends.
As I understand it I need the TLEN field from bam files (sam format field 9). However, all the alignment files from minimap2 have 0 in the TLEN field.
Does anyone know how to get pacbio alignments with the TLEN field? I am also happy using bwa or another aligner.
The command I used was minimap2 -ax map-pb {fasta} {reads} | samtools view -h -b - | samtools sort - > myfile.bam.

The TLEN for single-end reads is defined as 0 in the SAM specification, which is why this happens. You can still get the end position of the alignment by parsing the CIGAR string. An example using the pysam package in python is below:
import pysam

bam = pysam.AlignmentFile("myfile.bam")
for b in bam.fetch():
    print("This alignment spans {}:{}-{}".format(b.reference_name, b.reference_start+1, b.reference_end))

The reference_end value is automatically calculated from the CIGAR string.


text:
how to filter out cell after doing QC
I was wondering that the threshold that we set after doing QC in Seurat. there is an example
# GenePlot is typically used to visualize gene-gene relationships, but can be used for anything
# calculated by the object, i.e. columns in object@meta.data, PC scores etc.
# Since there is a rare subset of cells with an outlier level of high mitochondrial percentage
# and also low UMI content, we filter these as well
par(mfrow = c(1, 2))
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "percent.mito")
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "nGene")`

Then its filter cell command is under:
# We filter out cells that have unique gene counts over 2,500 or less than 200
# Note that low.thresholds and high.thresholds are used to define a 'gate'
# -Inf and Inf should be used if you don't want a lower or upper threshold.
pbmc <- FilterCells(object = pbmc, subset.names = c("nGene", "percent.mito"), low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))`

So I want to know whether the threshold setting is depending on the picture I show?
And the below is the Geneplot I have got base on my data. How can I set my thresholds?

Why have an upper bound on ngenes? Doesn't it indicate a cell is healthy if it is is able to express many genes?

Yes, you end up placing filters where things become inconsistent or similar to background. In the example, which appears to be from the Seurat tutorial, there's a downward plunge in nGenes vs. nUMI at around 200 genes and few believable cells with more than 2000-2500 genes, so the tutorials sets filters around there. The 0.05 mitochondrial contamination is a more arbitrary round number and is less supported by the plots (one would use more like 0.06 or 0.07 if one based this solely on the plots).

Thank you@Devon Ryan If my geneplot result shows below/Users/huali/Desktop/Screen Shot 2018-08-15 at 7.35.30 PM.png.                 How should I setting my threshold ?

What do you think would be good decent values given those plots? I don't want you to just blindly follow what I say :)

I just don't know how people set their thresholds, and want to figure out how they can do to do further analysis base on their QC. So I show the data that I got. Want to find something that can be dig out.

In your case percent.mito should be less than 0.1 and nGenes should be at least 1000.

From the graphs, what would be a good upper limit cut-off for nFeatures per UMI?


text:
Combine VCF files
I have a problem with using rbind to combine VCF files using the library VariantAnnotation from Bioconductor.
I am reading two VCF files, when I try to combine them in a certain order with rbind I'm getting an error. When I combine them in a different order, it works fine.
vcf.full <- rbind(vcf1, vcf2) ## Doesnt work

Error in as(from, to_class, strict = FALSE) :    no method or default
for coercing “CompressedCharacterList” to “DNAStringSetList”

vcf.full <- rbind(vcf2, vcf1) ## Work

header(vcf1)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

header(vcf2)
class: VCFHeader
samples(503): NA06984 NA06989 ... NA20828 NA20832
meta(2): META contig
fixed(2): FILTER ALT
info(27): CIEND CIPOS ... EX_TARGET MULTI_ALLELIC
geno(1): GT

dim(vcf1)
[1]  30 503

dim(vcf2)
[1] 149 503

The two VCF files seem identical. Any idea of why it happens? Or any suggestions on a different method to combine multiple VCF files using R?

Answer from @Kozolovska, converted from comment:
Found the problem, some of ALT were structural variant, so the list of vcf2@info$ALT was a charterer list and not a 'DNAStringSetList' after removing that variant it worked. Thanks.


text:
Visualize PLINK pedigree files
Is there a pedigree viewer that works with PLINK pedigree files? I saw an old post by Pierre, but I don't see the script on his blog.

We use R package kinship2 here: https://cran.r-project.org/web/packages/kinship2/index.html

We use kinship2 R package to visualize pedigree files (.fam plink files). For some reasons, some characters from plink ped aren't supported (like "-9" for unknown genre or affected status).
Here is the R script we use:
args<-commandArgs(TRUE)
pedPath<-args[1]
pdfFile<-args[2]

library(kinship2)
pedfile <- read.table(pedPath,sep="\t")
pedfile$V5 <- apply(pedfile,1,function(data){if(data[5] == 1 || data[5] ==2){data[5]}else{3}})
pedfile$affected <- apply(pedfile,1,function(data){if(as.numeric(data[6]) == 1){res <- data[6]}else{if(as.numeric(data[6]) == 2){res <- data[6]}else{res <- NA}}})

familyIDs <- as.character(unique(pedfile$V1))

pdf(file=pdfFile)
for(famID in familyIDs){
    ped <- pedigree(id=as.character(pedfile$V2[pedfile$V1==famID]),dadid=as.character(pedfile$V3[pedfile$V1==famID]),momid=as.character(pedfile$V4[pedfile$V1==famID]),sex=as.numeric(pedfile$V5[pedfile$V1==famID]),affected=as.numeric(pedfile$affected[pedfile$V1==famID]))
    plot.pedigree(ped,cex=0.7,align=TRUE)
    title(famID)
}

And we call it this way:
Rscript <kinshipScript.R> <plink ped file> <pdf output>

Now, if only I could convince people to include affected information in the actual .fam files. :-)

@MarkEbbert affected status is not really necessary but it still looks good on the pdf :)


text:
Stand alone chaining tool for existing blast hits?
I have a table of blast hits in a database that I'd like to chain. What is a good choice for a 'stand alone' chaining tool? Is it simple to implement in Perl, for example?
Example 'hit' data is here. Background is here

Do you know any existing stand alone chaining tools? How do you define good choice (fast, accurate, easy,..)? You have the description of the algorithm on page 19. And some implementations linked on page 38.

Ideally I'd like something I can understand in a short subroutine.

What is short? I find the question unclear to search for solutions or to build such short "subroutine".

I'm not able to implement the algorithm on page 19. I'm looking for an implementation of that algorithm in Perl or Python that I can feed my hits to.

it is easier to get help how to correctly implement an algorithm than get an algorithm implemented for you. If you could post your code and your problems in the question someone could help you to implement it. How wasn't the implementations on page 38 not helpful?


text:
How to map reads shorter than 32bp with minimap2?
I mapped fasta sequences to a reference. The fasta sequences range in length from 17bp up to several hundred bp. I used minimap2 with the following command: minimap2 -ax sr ref.fa seqs.fa, et cetera. Only the sequences greater than 32bp len mapped though.
samtools view -F 4 filename.sorted.bam | \
  cut -f10 | awk '{print length($1)}' | \
  sort | uniq -c | sort -k2 -n

   2304 32
  34931 33
   5443 34
   5343 35
   8811 36
   5485 37
   5417 38

Is it possible to make minimap2 map sequences shorter than 32bp? Do I need to use bwa aln or some other mapper?

Minimap2 is not the most optimal aligner in this case, it was designed for long reads, of which 32bp reads are clearly not. You could try lowering the  mnimizer k-mer length and with -k the default being 15, also try lowering the minimzer window size with -w.  However I suggest using bwa aln and bwa samse for all reads under 32bp and then merging the alignment later. Although I struggle to see what the utility of these short reads will be against a backdrop of longer reads?

They aren't really reads per se - just short sequences from another application that I need to align. Thanks for your comments.

Minimap2 doesn't work well with 32bp reads by design. Older mappers like bowtie1 and bwa-aln will be better.


text:
Do people today with a disease caused by single gene mutation that is inheritable - all share a common ancestor?
First question here so sorry if it is about both genetics and evolution.
So, if I get a disease which can be clearly identified on the gene level and all my descendants inherit it, after 500 years do all people with the disease share the common ancestor (given that it's only inheritable)?
I mean, my clear guess first would be "yes" but I really want to check it since I haven't read it anywhere.
Feel free to post materials on this too.
Thanks!

That can happen, but it's more likely that the disease will be present in more than one family. Diseases with a single, small genetic cause will be more likely to flare up elsewhere (assuming the disease isn't too bad, and it wouldn't be if it's been in the family for 500 years).
However, what you will find is that in your family that disease is almost always present together with some other inherited genetic feature that is nearby on the genome (i.e. an inherited haplotype block). This association is less likely to be shared with other families.


text:
Calculating read average length in a Fastq file with bioawk/awk
I found here this awk script:
BEGIN {
  headertype="";
  }
{
  if($0 ~ "^@") {
    countread++;
    headertype="@";
  }
  else if($0 ~ "^+") {
    headertype="+";
  }
  else if(headertype="@") { # This is a nuc sequence
    len=length($0);
    if (len>4) {
      readlength[len]++;
    }
  }
}
END {
  for (i in readlength){
    countstored+=readlength[i];
    lensum+=readlength[i]*i;
    print i, readlength[i];
  }
  print "reads read = "countread > "/dev/stderr";
  print "reads stored = "countstored > "/dev/stderr";
  print "read average length = "lensum/countstored > "/dev/stderr";
}

and I just wonder if it is possible to shorten it with bioawk?

It might be possible.  What problem do you have with having a long script? Why do you want a shorter script? Or do you want a faster script?

Did anyone get the correct script for finding out the average read length?

This can also be done with regular awk.
awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>

The NR%4==2 count the second line out of every block of 4.  length is a built-in that defaults to the length of the line, same as length($0).  In this case you can inject you custom printing to the END{} block but countread and countstored will always be the same since the averaging is done on the fly.

Note that, like OP’s code, this will only work on FASTQ files that don’t break the sequence and/or quality string across multiple lines.

That is correct.  However that practice is highly frowned upon so should not be an issue with most FASTQ files

I've yet to ever encounter a FASTQ file which has sequence or quality on multiple lines, I think for all intensive purposes assuming 4 lines per record should be fine.

This script is wrong because a quality string may start with @. With bioawk, it can be simplified to:
bioawk -c fastx '{ readlength[length($seq)]++; countread++ } END{...}'

The END{} block is the same as your original version.
EDIT: forgot -c fastx is the original answer. Thank @MattBashton for pointing this out.

Using the exact same END{} code from the OP fails in this example with bioawk: illegal field $(), name "seq"


text:
Efficiently aligning a lot of reads on the same small reference sequence
The context: I have a DNA-sequence coding for a protein, about 1500 bp in length. Using NGS, a lot of reads of (mutants of) this same sequence were acquired. All of these reads need to be aligned to the reference. We're talking about a lot of reads (100,000 - 1,000,000) of not too short length (120 - 300 bp). All of these reads can belong to different mutants of the template sequence, so the alignment is necessary to determine the exact sequence of every single mutant.
Currently, I'm just using Smith-Waterman-based local alignment to align every single read to the reference one by one. Yet I can't help but feel there might be a more computationally/time-efficient solution to this specific problem.
Maybe there exists an algorithm that's not very efficient for most alignment problems, but that becomes very worthwile if it has to map a ton of reads to the same place over and over again. For example, it might do some time-consuming operations on the short template that make it fast to align reads to it, but these operations aren't worth it if only aligning a couple of reads. That's just an idea, I don't know all the different techniques that are out there.
So, to recap: I want to align a lot (100,000 - 1,000,000) of 120-300 bp reads to the same short 1500 bp reference sequence. If anyone has any suggestions about an algorithm, or just a specific workflow, that is particularly suited to do this, it would be appreciated. I work in R so I can implement some things myself, it doesn't have to be a ready-to-use software package or anything like that. Thanks in advance!

Have you tried using some multi sequence alignments tools like T-coffe, or clustal? How many time are you willing to spend/how much accuracy do you want? Do you know which reads belong to each single mutant or is all in a pool?

@Llopis: yes, it is known which reads belong to the same mutant. I haven't used multi sequence alignment because I don't want to align all sequences to eachother, as I'm only interested in differences with the reference. As for time: at this moment it's still manageable to use common local alignment algorithms like Smith-Waterman, but if there's an approach that's just faster than this, it will help in the scale-up to even more high-throughput.

Have you tried using whole genome aligners like BWA and BOWTIE2 with your custom reference?  These are designed for NGS input and alignment of larger numbers or reads against a single reference.

@Bioathlete: No, I haven't tried it yet because BWA and BOWTIE2 were indeed designed for "whole genomes" while my reference is extremely small. I wondered if there were any methods that were commonly used in that case. Maybe not, and BWA or BOWTIE might still be the optimal choice, even for a small reference. That's as good an answer as any. I'll wait a little longer for extra suggestions, otherwise I'll consider the question answered. Thanks!

I have successfully used BWA and Bowtie2 for this sort of task, you may want to explore various parameters and how they affect the overall alignment rate.

@Bioathlete you may as well post that as an answer so the question isn't left hanging.

I would suggest using a whole genome aligner like BWA or BOWTIE2 and swapping out the whole genome reference with your specific reference.
These alignment algorithms are designed to align large numbers of NGS reads against a single reference.  They can be used with any reference including short ones like your example or small viruses.
To ensure that local alignment is correct, keep in mind that bowtie2 needs the --local while bwa should be run with on the BWA-MEM modality


text:
How to create a legend in FeaturePlot with do.hover at the same time?
you can see in my feature plot I can only see the cells from which sample but I cannot see the legend. And if I change the command I can see the legend but cannot see the sample identities.FeaturePlot(object = VarMerge, features.plot = c("Hcn4"), do.hover = T, no.legend = F,)
if I change the command to
FeaturePlot(object = VarMerge, features.plot = c("Hcn4"), no.legend = F

it will show like the second picture. So I just want to know whether in the first picture can I do something to the command and add the legend like in second picture?