text:
How do I add external labels and connecting lines to a Ternary Plot in R?
I have plotted a triangular plot,that three dimensions of which represent three herbal strategies. One dimension represents the strategy of C (competitive plant), the second dimension “S” (stress tolerant plants) and the third dimension”R” (ruderal plants), as you can see in the image below,
the points on it represent the plant species. Now I am trying to write name of  the plant species outside the triangle and connect with the arrow to the corresponding point. How do I do this work?
Many thanks in advance as it would help me a lot in my work.
TernaryPlot()
#Plot two stylised plots side by side, and plot data
par(mfrow=c(1, 1), mar=rep(0.3, 4))
TernaryPlot(atip='C%', btip='R%', ctip='S%',
            point='UP', lab.cex=0.8, grid.minor.lines=0,
            grid.lty='solid', col='#FFFFFF', grid.col='GREY',
            axis.col=rgb(0.1, 0.1, 0.1), ticks.col=rgb(0.1, 0.1, 0.1),
            padding=0.08)
data_points <- list("Bromus dantonia" = c(47,   59, 149),
                    "Calamagrosis psoudo phragmatis" = c(90,    102,    63),
                    "Carex diluta" = c(109, 64, 82),
                    "Carex divisa" = c(96,  99, 59),
                    "Carex pseudocyperus" = c(130,  71, 54),
                    "Carex stenophylla" = c(97, 98, 59),
                    "Catabrosa aquatica" = c(100,   5,  150),
                    "Centaurea iberica" = c(124,    85, 46),
                    "Cirsium hygrophilum" = c(158,  42, 55),
                    "Cladium mariscus" = c(159, 96, 0),
                    "cod2" = c(54,  82, 119),
                    "Cynodon dactylon" = c(121, 54, 80),
                    "Eleocharis palustri" = c(124,  100,    31),
                    "Epilobium parviflorum" = c(67, 80, 107),
                    "Eromopoa persica" = c(83,  15, 157),
                    "Funaria cf.microstoma" = c(8,  0,  247),
                    "Glaux maritime" = c(4, 196,    55),
                    "Hordeum brevisubulatum" = c(76,    70, 109),
                    "Hordeum glaucum" = c(40,   79, 136),
                    "Inula britannica" = c(95,  108,    51),
                    "Juncus articulatus" = c(107,   79, 69),
                    "Blysmus compressus" = c(81,    127,    47),
                    "Juncusinflexus"= c(149,    106,    0),
                    "Medicago polymorpha" = c(60,   86, 109),
                    "Mentha spicata" = c(150,   23, 82),
                    "Ononis spinosa" = c(66,    112,    77),
                    "Phragmites australis" = c(234, 0,  21),
                    "Plantago amplexicaulis" = c(108,   83, 64),
                    "Poa trivialis" = c(90, 28, 138),
                    "Polygonum paronychioides" = c(20,  12, 223),
                    "Potentila reptans" = c(106,    41, 108),
                    "Potentilla anserina" = c(105,  58, 91),
                    "Ranunculus grandiflorus" = c(129,  25, 101),
                    "Schoenus nigricans" = c(143,   91, 21),
                    "Setaria viridis" = c(10,   7,  238),
                    "Sonchus oleraceus" = c(178,    0,  77),
                    "Taraxacum officinale" = c(117, 28, 110),
                    "Trifolium repens" = c(94,  4,  157),
                    "Triglochin martima" = c(63,    96, 95),
                    "Veronica anagallis-aquatica" = c(55,   37, 163)
)
AddToTernary(points, data_points, pch=21, cex=1.2,
             bg=vapply(data_points,
                       function (x) rgb(x[1], x[2], x[3], 128, maxColorValue=255),
                       character(1))
)

AddToTernary(text, data_points, names(data_points), cex=0.8, font=1)

Hmmm I'm not familiar with AddtoTernary function, I can see that you've placed a function within the "AddtoTernary" function. The traditional ggplot2 way of doing this would simply be controlled by the theme
ggtern(data=x,aes(x1,x2,x3)) +
  geom_point(fill="red",shape=21,size=3) +
  theme_tern_bw()

This will give an arrow for each axis, but additional code is required for different colours (I'd have to look up my past work).
I think what you might want is a colour key? In this case ggtern will have a "key" function because ggplot2 has one.Here, your colouring is based on species not on a colour gradient as to how far to the corner of a graph you are.
Alternatively you can easily put names besides each point in ggplot2, but again this is traditional style code.
Also you can add different shapes to the graph, which might be useful instead of colours.
Arrows and names together peronsally I've not seen on a ggplot2 graphic and I'd have not idea how to do it.
FYI I don't think we need all the data just a few points fine, although it shows why you prefer arrows - personally I would use colour coding. In "traditional" ggplot2 coding the graphing is built up in layers, so if additional functionality is required it is simply appended to the exisiting command. Thus it is quite easy to make modifications.

Hi, Thanks for answering Unfortunately, I was not able to draw the chart with the ggplot2 package.I want to write the species name outside the triangle as shown below and connect it to the points inside the triangle with the arrow.how do I draw this ternaryplot? Thank you, I really appreciate your help.


text:
What is "aligned sequences" and "consensus sequence" in the context of sequence logo? How to compute these?
In bioinformatics, a sequence logo is a graphical representation of the sequence conservation of nucleotides (in a strand of DNA/RNA) or amino acids (in protein sequences).
so, What is "aligned sequences" and "consensus sequence"?
how can I get these sequences based on following data?

Figure 2.5 from the book "Machine Learning - A Probabilistic Perspective"
more info about this figure at here

Is this homework?

@KamilSJaron from the book "Machine Learning - A Probabilistic Perspective"

@KamilSJaron https://stats.stackexchange.com/q/419575/247356

@KamilSJaron what does "AT" in "be made of AT" mean?

Panel (a) shows a sequence gap-free alignment. Each row corresponds to a contiguous 15 base pair sequence of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.). Each column corresponds to a particular position/residue (e.g. column 1 is the 1st position of the alignment etc.).
Panel (b) shows the sequence logo. As noted in the comments, this does not correspond to the alignment in (a). Instead, here is the corresponding logo:

From this you can intuitively understand that at the 1st sequence position you have only As and Ts and their frequency corresponds to the y-axis. For the 3rd position, you have only As, so only A is represented here in the logo etc.
The consensus sequence is the most frequent letter at each position in the alignment, i.e. [AT]TAGC etc. where [AT] means A or T since they are jointly high in frequency (5 each).

source: Stack Exchange

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

sort: votes


text:
How to use -csvStats option?
I'm trying to get CSV files from snpEff -csvStats option and got stuck. I tried the code below and only got the vcf file with no contents. How can I use -csvStats option properly?
snpEff -Xmx4g -csvStats -v hg19 08_FILTER/NA12878.Filtered.Variants.vcf > test.ann.vcf

    00:00:00  SnpEff version SnpEff 4.3t (build 2017-11-24 10:18), by Pablo Cingolani
    00:00:00  Command: 'ann'
    00:00:00  Reading configuration file 'snpEff.config'. Genome: '08_FILTER/NA12878.Filtered.Variants.vcf'
    00:00:00  Reading config file: /home/ubuntu/01_NA12878/snpEff.config
    00:00:00  Reading config file: /home/ubuntu/miniconda2/share/snpeff-4.3.1t-2/snpEff.config
    java.lang.RuntimeException: Property: '08_FILTER/NA12878.Filtered.Variants.vcf.genome' not found
      at org.snpeff.interval.Genome.<init>(Genome.java:106)
      at org.snpeff.snpEffect.Config.readGenomeConfig(Config.java:681)
      at org.snpeff.snpEffect.Config.readConfig(Config.java:649)
      at org.snpeff.snpEffect.Config.init(Config.java:480)
      at org.snpeff.snpEffect.Config.<init>(Config.java:117)
      at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
      at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
      at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
      at org.snpeff.SnpEff.run(SnpEff.java:1183)
      at org.snpeff.SnpEff.main(SnpEff.java:162)
    00:00:00  Logging

00:00:01  Done.

 snpEff -Xmx4g -v hg19 -csvStats 08_FILTER/NA12878.Filtered.Variants.vcf > tes2.ann.vcf

00:00:25  Genome stats :
#-----------------------------------------------
# Genome name                : 'Homo_sapiens (USCS)'
# Genome version             : 'hg19'
# Genome ID                  : 'hg19[0]'
#     'M' 16571   Vertebrate_Mitochondrial
#     'Un_gl000226'   15008   Standard
#     '18_gl000207_random'    4262    Standard
#     'MT'    1   Vertebrate_Mitochondrial
#-----------------------------------------------

00:00:27  Predicting variants

Here is the snpEff.config file. I installed the snpEff at /miniconda2/share
.
# Databases are stored here
# E.g.: Information for 'hg19' is stored in data.dir/hg19/
#
# You can use tilde ('~') as first character to refer to your home directory.
# Also, a non-absolute path will be relative to config's file dir
#
#---
data.dir = ./data/

#---
# Database repository: A URL to the server where you can download databases (command: 'snpEff download dbName')
#---
database.repository = http://downloads.sourceforge.net/project/snpeff/databases

#---
# Latest version numbers. Check here if there is an update.
#---
versions.url = http://snpeff.sourceforge.net/versions.txt

#-------------------------------------------------------------------------------
# Third party databases
#-------------------------------------------------------------------------------

# Databases for human GRCh38 (hg38)

database.clinvar.GRCh38                 = ./db/GRCh38/clinvar/clinvar-latest.vcf.gz
database.repository.clinvar             = ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar-latest.vcf.gz

database.dbsnp.GRCh38                   = ./db/GRCh38/dbSnp/dbSnp.vcf.gz
database.repository.dbsnp.GRCH38        = ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz

database.dbnsfp.GRCH38                  = ./db/GRCh38/dbNSFP/dbNSFP.txt.gz

# Databases for human GRCh37 (hg19)
database.clinvar.GRCh37                 = ./db/GRCh37/clinvar/clinvar-latest.vcf.gz
database.repository.clinvar.GRCh37      = ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar-latest.vcf.gz

database.dbsnp.GRCh37                   = ./db/GRCh37/dbSnp/
database.repository.dbsnp.GRCh37        = ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b141_GRCh37p13/VCF/

database.gwascatalog.GRCh37             = ./db/GRCh37/gwasCatalog/gwascatalog.txt
database.repository.gwascatalog.GRCh37  = http://www.genome.gov/admin/gwascatalog.txt

database.dbnsfp.GRCh37                  = ./db/GRCh37/dbNSFP/dbNSFP.txt.gz

#-------------------------------------------------------------------------------
# Loss of function (LOF)
#-------------------------------------------------------------------------------

# It is assumed that even with a protein coding change at the
# last 5% of the protein, the protein could still be functional.
lof.ignoreProteinCodingAfter  : 0.95

#  It is assumed that even with a protein coding change at the
#  first 5% of the protein:
#       "..suggesting some disrupted transcripts are
#       rescued by transcriptional reinitiation at an
#       alternative start codon."
lof.ignoreProteinCodingBefore : 0.05

# Larger deletions removing either the first exon or more than
# 50% of the protein-coding sequence of the affected transcript
lof.deleteProteinCodingBases : 0.50

#-------------------------------------------------------------------------------
# Codon tables
#
# Format:   It's a comma separated "codon/aminoAcid[+*]" list
#           Where 'codon' is in uppper case, aminoAcid is a one letter
#           code, '+' denotes start codon and '*' denotes stop codon.
#
# References:   http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi
#               ftp://ftp.ncbi.nih.gov/entrez/misc/data/gc.prt
#-------------------------------------------------------------------------------

Hmmmm. The message Genome: '08_FILTER/NA12878.Filtered.Variants.vcf' is confusing me. That's not in this snpEff.config file. Is it in the other config file? The terminal output includes /home/ubuntu/01_NA12878/snpEff.config and /home/ubuntu/miniconda2/share/snpeff-4.3.1t-2/snpEff.config. If the VCF file isn't in either config file, my next guess would be that the command-line arguments are declared incorrectly. I don't know enough about snpEff to troubleshoot that without more context, though. :(

@DanielStandage I've got the one more issue with the installation, I thought it might be something with snpEff installed in miniconda2, so I deleted the snpEff and re-installed it in the prior directory of it, I realized it's not a good idea so I removed the snpEff again and installed the snpEff at /share/miniconda2 again. And now I'm getting error snpEff: command not found, how can I make this working?

@DanielStandage I mean, how can I make just 'snpEff' code working without java jar? It worked before I re-installed it but now I think something is twisted. I don't have an idea cause    snpEff -Xmx4g -v hg19 08_FILTER/${sample}.Filtered.Variants.vcf > 09_ANN/${sample}.ann.vcf
worked without -csvStats option well.

Had the same issue, but found that I need to provide a filename to use the csvstats option
snpEff -Xmx4g -csvStats stats.csv -v hg19 08_FILTER/NA12878.Filtered.Variants.vcf > test.ann.vcf

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


text:
Why do ten rows (Figure_1) correspond to 2 bits (Figure_2) in a sequence logo?
Following this question, I'm confused with the computation of sequence logo
Following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure_1)"

here is the corresponding sequence logo (Figure_2).

ten rows represent sequences of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.)
Each column corresponds to a particular position.
Why do ten rows (Figure_1) correspond to 2 bits (Figure_2)?

You can find on e.g., wiki how a sequence logo is created.

Bits are not frequencies. If a position only contains an A (position 3 for example) then you would need 2 questions (bits) to derive that information without a priori knowledge.

Is it a G or C? if no > then it is a A or T
Is it T? If no then it is an A

Position 1 can be derived by 1 question only:

is it a G or C?  If no then it is an A or T

In this case the number of bits can be calculated using the frequenci)es.
For 4 bases, A only (100%):
4 * -0.25 * log2(0.25) - (-1 * log2(1)) = 2 (bits)
For 8 letters, A only (100%):
8 * -1/8 * log(1/8)- (-1 * log2(1)) = 3 (bits)
Mixed base possibility (A/T 50% each) with 4 bases:
4 * -0.25 * log2(0.25) - (-0.5 * log2(0.5)+-0.5 * log2(0.5)) = 1 (bit)
As mentioned in the comments, you can follow the wiki link for an example.
Or, an example with a little more explanation from the biology SE site.
sequence logo calcuation
Hopes this helps

Thanks for your answer. In this specific case, 4 bases is being used, since there are a,t,g,c 4 Nucleotides in total, right?

correct. That is this part: 4 * -0.25 * log2(0.25) of the equation.


text:
How to retrieve the subcellular location info from uniprot?
I have a list of proteins with their gene names/symbols. I would like to identify which of them are extracellular proteins or plasma membrane proteins.
I can get this information in UniprotKB by searching the proteins one by one. Because I have hundreds of proteins, I wonder whether I can do a systematic search of the UniprotKB in RStudio. I tried the UniProt.ws select method. But I couldn't even get the example code in the package work. I downloaded the uniprot_sprot.dat file. But I don't know how to read this file correctly in RStudio. I used tab-delimited format. It turned out there is only one variable in the file.
Can anybody tell me what I should do to solve the problem? I am pretty new to any computer languages.

Can you give us an example ID and where on the UniProt website you see the field you're looking for? Examples always help make the question more specific and grounded.

Welcome to the site. I edited a bit the question to make it easier to find it and read. You could use gene ontology cellular component to know where a protein can end. Also, show how did you try to read the file and why did you download that file? I couldn't find anywhere that it has the information you seek.

Since you are already on UniProtKB, the simplest way will be if you add this information to the output table and download it as a that you can import and filter in R.
To get additional annotations, click the button labelled "Columns" on top of the UniPortKB website.

This allows you to add a wide variety of annotations such as subcellular location to the output. Tick the checkmark for subcellular location and save on the top right.

Then you can download the proteome-wide annotation for example as a tab-separated list using the "Download" dialog next to the "Columns" selection. This list can be imported into R.

Hi PPK,Thanks for the information on how to download the subcellular location info from UniProtKB. The file seems very large to be downloaded. I am in the process of downloading.

After my post of the question, I went to try the UniProt.ws package in the RStudio again. I realized that there is a mistake in the example code in the package. The up object is from a human database, but the keys the example provides are the ENTRZID from mouse proteins. When I tried human ENTRZID numbers, it worked. The only problem for the select method in the UniProt.ws is that it can't retrieve the protein information for a large batch of proteins. The maximum number that worked for me is about 80 proteins per time.


text:
How to compute the Shannon entropy for a strand of DNA?
I'm confused by the computation of sequence logo. Wikipedia gives a process about this without a concrete example.
Let's just consider DNA, so there are 4 bases (nucleic acids).
The following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure 1)"

The Shannon entropy of position $i$ is:
${\displaystyle H_{i}=-\sum _{b=a}^{t}f_{b,i}\times \log _{2}f_{b,i}}$
Where ${\displaystyle f_{b,i}}$ is the relative frequency of base
This post is computing position 3, where it seems that the relative frequency of a is 100%, and the relative frequency of 3 other bases are 0%.
Now, how do you calculate the Shannon entropy of position 3? It seems to be 0, which is apparently incorrect, where am I wrong?

Note that the frequency is always in the interval [0,1] instead of [0,100]. H is 0 for columns with only one symbol.

Why do you think the entropy of 0 is incorrect? It intuitively makes sense, as there is no uncertainty about the base at position 3, and thus there is no entropy.
However, what is plotted in a sequence logo isn't the entropy, but rather a measure of the "decrease in uncertainty" as the sequence is aligned. This is calculated by taking the entropy at this position if we randomly aligned sequences ($H_g(L)$), and subtracting from it the entropy of the alignment ($H_g(s)$):
$$
R_{sequence}=H_g(L) - H_s(L)
$$
The entropy at position 3 based on a random alignment is calculated by assuming there are 4 equally likely events (one for each base), and thus:
$$
\begin{align*}
H_g(L) & = -((1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2)) \\
H_g(L) & = 2
\end{align*}
$$
Notably, $H_g(L)$ will always be 2 when dealing with nucleotide sequences.
So in your example, we have:
$$
\begin{align*}
R_{sequence}&=H_g(L) - H_s(L) \\
R_{sequence}&=2 - 0 \\
R_{sequence}&=2 \text{ bits of entropy}
\end{align*}
$$
Finally, to work out the height of each base at each position in the logo, we multiply the frequency of that base by the overall information gain $R_{sequence}$. In this case the frequency of the base A at position 3 is of course 1:
$$
\begin {align*}
H(b,l) &=f(b, l) \times R_{sequence} \\
H(A, 3) &= f(a, 3) \times 2 \\
&= 1 \times 2 \\
&= 2
\end{align*}
$$
Source: Sequence logos: a new way to display consensus sequences.
 (Thomas D. Schneider, R. Michael Stephens, 1990)
Note: I've excluded the correction factor $e(n)$ from the overall information-gain calculation for the sake of simplicity. Refer to the paper above for an explanation of this.

Thanks for your answer. What does "g" mean/represent/come from in $H_g(L)$

I'm not sure. My guess is g=genome and s=sequence(ing) but that's pure conjecture since the paper doesn't spell it out

Thanks for your reply. And does L mean/represent/come from location since the paper says "at various positions (L)".

Yes, that would be my assumption. This is important to include in the equation because it's only calculating the information at a single position. The authors propose another equation for an entire sequence (which is just summing these per-position scores)

Thanks a lot. And last question, does R come from/represent/mean acids?

I don't know. But also I think you're putting too much importance on the authors' choice of variables. You could write any of these equations with different letters and they would be just as correct.


text:
Passing around complex metadata in Nextflow
I have a number of samples (from the 1000 genomes) that I want to analyse using Nextflow. However, I need to keep track of their metadata, fields like sample ID, analysis type etc. So I have a basic metadata object like this for each sample that I then pipe into a channel:
            metadata = [
                sample: split[0],
                mapped: split[1],
                platform : split[2],
                alignment: split[3],
                tsi: split[4],
                type: split[5],
                date: split[6],
                filetype: split[7..-1],
                url: file(bam)
            ]

Then, in the first process, I attempt to use this metadata object:
process bamToFastq {
    input:
        val bam

    container 'quay.io/biocontainers/samtools:1.8--4'
    cpus 4
    disk '50 GB'

    output:
        file('r*.fastq') into fastq

    """
    samtools fastq -0 r0.fastq -1 r1.fastq -2 r2.fastq -@ ${task.cpus} $bam.url
    """
}

However, I already run into some issues:

Nextflow doesn't realise that metadata.url is a file, since I declared the whole object using val, so I'm guessing it won't be correctly localized (based on the documentation).
I want pass down a new metadata object for each of my generated FASTQs, that's the same as the input one, but I want to add in the FASTQ as a URL, and the read number of each FASTQ. Where/how can I do this? Nextflow doesn't provide any section of the process for structuring the output data. If this were Python I would do the following:

    metadata = bam.copy()
    metadata['fastq'] = file
    metadata['read_number'] = file.name[1]

Is there a good way to pass around structured metadata like this in nextflow? Will this solve either of my problems?


text:
Subset a multisample VCF file
I need some help with GATK. I have a multisample vcf file containing 18125 samples. After performing PCA I have around 303 samples as outliers and I wanted to remove these samples from my multisample vcf file.
I tried doing it using the GATK SelectVariant tool. It does some stuff and generates a VCF file which is a bit smaller in size than the original file, but in the end the total count of samples remains the same. As I have no prior knowledge of GATK (working with it for the first time) I am just wondering if you guys can help me out in fixing this issue.
The command that I used was
gatk SelectVariants -V filtered.vcf -xl-sn remove.txt -O test.vcf

remove.txt contains the sample names in each row (1 sample per row).
After running this command it generated a test.vcf file, for which I counted samples using:
bcftools query -l test.vcf | wc -l

However, I got exactly same number of samples as in the original VCF file filtered.vcf.

Based on the way you've called the command (using -xl-sn), I assume you are using GATK 4 and not GATK 3.
If this is the case, the documentation for the -xl-sn flag is available here: https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_variantutils_SelectVariants.php#--exclude-sample-name.
By default, the -xl-sn SAMPLE means "exclude the sample called SAMPLE from my VCF". In your case, you're providing an entire list of samples, in a file, which is not the default behaviour. In order to tell GATK that you're providing a file and not a sample name, you need to ensure the file you provide has the extension .args.
From the documentation:

File names must use the extension ".args"

ok cool it worked thanks mate (Y)


text:
How to demultiplex a mix of single-indexed and dual-indexed samples
The problem
If I have a sample sheet that contains both single-indexed and dual-indexed samples, I can split it up into two sample sheets and then run bcl2fastq on each one. However, when doing this, large Undetermined fastq files are generated. E.g., when processing the single-indexed samples, all the dual-indexed samples go to Undetermined. And when processing the dual-indexed samples, the single-indexed samples go to Undetermined.
Additionally, because they are being processed separately, if any single index A is part of a dual index A+B, then when processing the single indexes, an A+B may be mistaken for an A, so it would seem that they need to processed simultaneously to avoid this mis-assignment.
The question
Given a sample sheet and directories of BCL files, how can such a set of sequencing data be demultiplexed correctly, either using bcl2fastq or the Picard tools?
To put it another way, I want to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. It can be assumed that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.
Attempts at a solution
Using bcl2fastq directly
If a sample sheet of a mix of indexes is given to bcl2fastq (v2.17) directly, it produces the error
ERROR: bcl2fastq::common::Exception: Success (0): .../bcl2fastq2/src/cxx/lib/layout/BarcodeCollisionDetector.cpp(127):
Throw in function void bcl2fastq::layout::BarcodeCollisionDetector::validateNewBarcodeSizesAgainstExisting(const std::vector<long unsigned int>&) const
Dynamic exception type: boost::exception_detail::clone_impl<bcl2fastq::layout::BarcodeCollisionError>
std::exception::what: Barcodes have an unequal number of components.

Barcodes have an unequal number of components.

Using picard
It seems that it should be possible using the Picard tools, but I have not found a way to set up the inputs to ExtractIlluminaBarcodes and IlluminaBasecallsToFastq that processes this configuration correctly.
The syntax for supporting multiple index configurations is not entirely clear. But when using various combinations of N and '*' on multiplex_params.tsv and barcodes.txt required by the Picard tools, a large Undetermined file is still produced and actual sample fastq files are tiny, indicating it is not processing them correctly.
Using multiple calls to bcl2fastq
As indicated in the discussion above, this suffers from the problem of having all the dual indexed samples going to the "Undetermined" file when processing the single-indexed samples, or vice versa, and creates two output directories of reports and stats, which must be merged.
Padding empty indexes with N
By padding, I mean converting ACGTACGT to ACGTACGT+NNNNNNNN so that the single-index samples in the sample sheet are "dual" as well. This strategy seems that it would work except there is a bug in bcl2fastq that it treats "N" literally instead of as a wildcard. See the release notes for details.

Welcome to the site. Sorry but I don't understand which files do you have. Could you clarify what are your files? Or do you have a mix of single-indexed and dual-indexed and you can't differentiate between them?

@llrs Hi, thanks for the question. In this case, I'm trying to demultiplex a single sequencing run that contains both single-indexed and dual-indexed samples. You can assume that the indexes are sufficiently distinct such that any sample's index configuration is different from any other. But assuming that the different index configurations are not segregated to particular lanes of the sequencer, the question is how to demultiplex the files correctly such that both the single-indexed and dual-indexed samples are recognized.

Ok, now I understood the question, but I'm sorry I can't help with it. Perhaps you could [edit] your last comment in the question for other users.

I'm late to the party but simply run cellranger mkfastq twice, with different arguments for the mask and with --filter-dual-index for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace consuming, but it can be done. I would use some Unix tools for that but there are probably smarter ways.
EDIT: I mistakenly said bcl2fastq, this must be cellranger mkfastq. (thanks AmadeusDrZaius)

Thanks for your suggestion. I only see that "--filter-dual-index" option available for cellranger mkfastq, rather than bcl2fastq itself. Is that where you saw the option?

ah yes, that is true, sorry. I guess you would have to do an additional filtering step on the FASTQ file, the second word on the ID lines is 'index1+index2' for properly demultiplexed reads.

Have you just tried giving bcl2fastq one sample sheet with a mix of single and dual indices?  I don't think what you are trying to do is a problem.

if any single index A is part of a dual index A+B,

Well, it might be a problem if you did that.  That was poor planning.

That was the first thing tried but unfortunately bcl2fastq will give the error "Barcodes have an unequal number of components" when you try to do that (using v2.17). I'll add a note to my question.

Are you sure that the single index samples don't have a site for the second index to sit down and read?  My lab does some samples with single indicies and some with double, but they always prep with libraries with kits meant for dual indexing, so the second index is always there, even if I ignore it.

That's a good idea that might be a good workaround in some cases, but I think the general answer is that they are not always physically dual-indexed. The use of "N" as described in the bcl2fastq documentation would work if not for the fact that a bug in bcl2fastq causes it to treat "N" as a literal base "N" instead of as a wildcard. Given that, I'm guessing that the Picard tools are more likely to have a full solution for this particular situation.

I don't think it's going to be trivial for software to figure out the difference between "second index doesn't match what I expect" and "second index isn't really here".  You might have to have bcl2fastq export the indices as fastqs, so you can use the quality score to parse yourself.

That's an interesting point. I expected that it would just make a "best guess", but I can see why it wouldn't. It sounds like that would primarily only be an issue in the case of "A" and "A+B" that I described above. If the case of "A" and "A+B" is avoided, it should be able to tell the difference. Thanks for your thoughts.

I would be surprised if the Experiment Manager software would allow a run configuration that was un-processable.  I take it the instrument was run without a sample sheet made by Experiment Manager?


text:
what does a sequence look like before alignment?
I am confused with Sequence alignment, which is

a way of arranging the sequences of DNA, RNA, or protein to
  identify regions of similarity that may be a consequence of
  functional, structural, or evolutionary relationships between the
  sequences.

to understand this easily, I would like to use this image (image_1)

the image above (image_1) is the primary part of nucleic acid structure (primary, secondary, tertiary, and quaternary), shown below (image_2)

is the sequence gcgc... an aligned sequence? if yes, what does the sequence look like before alignment?

Sequence alignment typically refers to the primary sequence. That is, we ignore the complex 3D structure of the molecule and its spatial arrangement in the cell. Instead, we use a symbolic representation composed of symbols like A, C, G, T, and U. Alignment is a conceptual tool that, as you pointed out, is used to find regions of similarity between sequences. Therefore, it only makes sense in the context of two or more similar sequences.
By way of example, consider a simple case of two sequences X and Y. This is what they would look like before alignment. The similarities between the two sequences are unclear.
X: ATTGACATAGGCACTTCGACTA
Y: GACATCTTCGGCTAA

This is what they would look like after alignment. The similarities (and differences) are now much easier to see.
X: ATTGACATAGGCACTTCCACTA-
      |||||     ||||  |||
Y: ---GACAT-----CTTCGGCTAA

This arrangement identifies regions of similarity that may be related through evolution. This doesn't always mean perfect agreement though. As this example shows, sometimes mismatches or gaps are allowed so that long stretches of identical matches are lined up correctly. Sometimes, for very closely related sequences, the optimal alignment doesn't require any gaps at all. This seems to be the case with some of your recent questions about sequence logos.
Note: Sequence alignment is a fundamental concept in bioinformatics, and there is a tremendous amount of research and development on the topic. It's hard to imagine doing any amount of serious DNA sequence analysis responsibly without a basic understanding of sequence alignment. I would recommend searching Google and/or Google Scholar for intro-level reviews of the topic.

There are a lot of really good materials by Ben Langmead (the guy behind Bowtie  and HISAT) on sequence alignment, what it means, and how it's done. I'd also recommend most things written by Heng Li, if you're looking for places to start learning about sequence alignment


text:
Intersecting multiple data frames
I have to intersect multiple data frame to find common elements. For visualization purpose i can use upset Plot but it doesn;t return a object containing the intersection object .I tried jvenn but at max it can take six sets to intersect .
Im able to do the interesction but i would like to know if i can do it more efficiently ,which apparently not sure how to do this is my code as of now where i read each file as data frame, covert them to list and pass them to do intersection .
NewList <- split(GO0000122, f = seq(nrow(GO0000122)))
NewList1 <- split(GO0006351,f=seq(nrow(GO0006351)))
NewList2 <- split(GO0006355,f=seq(nrow(GO0006355)))
NewList3 <- split(GO0006357,f=seq(nrow(GO0006357)))
NewList4 <- split(GO0006366,f=seq(nrow(GO0006366)))
NewList5 <- split(GO0030154,f=seq(nrow(GO0030154)))
NewList6 <- split(GO0045892,f=seq(nrow(GO0045892)))
NewList8 <- split(GO0045893,f=seq(nrow(GO0045893)))
NewList7 <- split(GO0045944,f=seq(nrow(GO0045944)))

a <- intersect_all(NewList,NewList1,NewList2,NewList3,NewList4,NewList5,NewList6,NewList7,NewList8)

b <- do.call(rbind.data.frame, a)

Can i read the files and pass them into list and do the same that would help me from doing lot of things manually
I can read the files from the directory where i have my files but not sure how to pass them to make it as list.
filenames <- gsub("\\.csv$","", list.files(pattern="\\.csv$"))
for(i in filenames){
  assign(i, read.csv(paste(i, ".csv", sep="")))
}

The object filenames is character here .It just returns me list of files.
So how to make my question bit more clear , i would like to read the list of file from the directory and then pass these to perform intersection. Any suggestion or help would be really appreciated .

I've been cleaning some of your post for punctuation, but this one I am not sure how to change it. Could you please [edit] and correct the comas and i to I and other punctuation errors? Also, if you could clarify what do you mean by intersection here I think that your problem might be easier to solve if you explain what is your goal with this (there might be an out of the box solution)

oh sorry typing partly from phone adds to the problem. My goal is to get set of common elements from multiple data frame where their is only single column consist of genes.So as of now i read each file,then pass those as list, then i do intersect and then turn it into dataframe. So now each time i have to read files if i have 12 files i have to read 12 times .Is there a way to put it in a loop and do the same what i have done above.

I think it would be it easier to work with character vectors than lists (it will at least make it easier to see them). But what have to do the files with the intersection, if any first you read and then you intersect, not the other way around or I am missing something here?

it would be easier if i add my files which im using  https://drive.google.com/open?id=1laieKec1sr4rWa4SIhUvHcLSerrWf4OC

Easier for you but not for me. I don't understand what have to do this files with your problem

okay may be its too confusing but can you let me know if i can how can i read all the csv files in a loop and make them into data-frame that would be easier for me to do the rest

This question is out of the scope of this site. You can find your solution at StackOverflow, specifically this question

sure this helps

Change dataframe to load all the files into a list:
wd <- "C:/all_your_dataframes_in_one_folder"
setwd(wd)
Then you get the list of files you want from your folder
files <- list.files()
Put into a list with a loop
dflist <- list()
for(i in 1:length(files)){
 dflist[[i]] <- read.csv(files[i])
}

Then use tidyverse solution
library(tidyverse)
map(dflist, ~.$GeneID) %>% #creates a list just of the column of interest
    reduce(intersect)  #applies the intersect function cumulatively to the list

If you need a package it would be better if you state which ones. Also with readr package you can read already a bunch of csv files in one single call

I mentioned tidyverse, but I just added library(tidyverse) in the example to be extra clear. Also, I felt as though a loop example would be better for explaining how to construct them..

" a loop example would be better for explaining how to construct them" this is what im looking for even though i dont use much loop etc because Im sure or my concepts are pretty unclear about it.

@h3ab74  i ran your code    map(dflist, ~.$GeneID) %>% reduce(intersect)
I get this NULL .Am i doing something wrong

the part of the code with ~.$GeneID is supposed to be the column name of all your dataframes in the list. So, if the column with the gene names in your dataframes does not have the name 'GeneID', change the GeneID part in the code to one that works with your dataframes

In my case the column name of all my data frame header is Gene, so i ran this
map(dflist, ~.$Gene) %>%
reduce(intersect)
character(0)  I get character(0)

Can you edit your question to include an example of how the dataframes are structured so I can see why you're getting nothing? Also, please note that this solution will get you the genes that are similar in all of the dataframes

@h3ab74 these are my files https://drive.google.com/open?id=1laieKec1sr4rWa4SIhUvHcLSerrWf4OC my data frame which i import is single column with a header name Gene

it worked ,it had name issues.


text:
BLAST server download
I discovered that Phytozome's BLAST results look quite nice (please see below).
By any chance, does anyone know where to download this BLAST server in order to install it locally for organisms which are not part of Phytozome?

Thank you in advance,

I don't know about the BLAST server you're referring to unfortunately.
But another BLAST server, SequenceServer looks quite nice too. And has a bunch of helpful visualisations that you can download directly in png or svg. Source code and installation instructions for this one are at http://sequenceserver.com
My opinion is biased because my lab makes SequenceServer. But there is also social proof: dozens and dozens of community websites use it, and it get cited regularly.
Apologies for spamming a bunch of screenshots. I think they're helpful to understand what SequenceServer can do & I'm proud of what we've built.
(it's open source so you can install it on your own server, or use our cloud BLAST hosting service).
Overview of how hits align to query

Lengths of queries and hits

Kablammo pairwise comparison of hit-query alignment

Circos overview of BLAST results

Search interface (automatically detects query type & chooses BLAST algorithm)

Example BLAST result overview


text:
Multiple line chart
I am trying to create a multiple coloured line graph based on six lines.
Essentially I have data that identifies different units of vegetation by the variable” NDVI (vegetation index)” I want to create a plot with differently coloured lines for each units.
The following is a screen shot of the data structure and my code.
It is plotting all of the data. But the time data on X-axis is disrupted and not ordered. What am I missing? And also how do I get legend for my chart?
library(ggplot2)
library(tidyverse)

df2 <- data.frame(Time<-c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  A<-c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  B<-c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  C<-c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  D<-c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  E<-c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  F<-c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))
head(df2)

ggplot(df2, aes(x = Time, group = 1)) +
  geom_line(aes(y = A), colour="blue") +
  geom_line(aes(y = B), colour = "red") +
  geom_line(aes(y = C), colour = "green") +
  geom_line(aes(y = D), colour = "black") +
  geom_line(aes(y = E), colour = "yellow") +
  geom_line(aes(y = F), colour = "brown") +
  ylab(label="NDVI") + xlab("Time")

Couple of things regarding your code:

It would be good to have your data frame in long format for ggplot2

For having the x-axis ordered chronologically, you have to specify it manually, otherwise it is ordered alphabetically.

Then, here is my solution:
library(ggplot2)
library(reshape2)

# I have substituted your <- by = and changed the name of the groups
df2 <- data.frame(Time = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"),
                  g_A = c(0.209,0.225,0.237,0.312,0.377,0.468,0.499,0.492,0.45,0.43,0.405,0.405,0.378,0.358,0.337,0.297,0.179,0.305,0.241,0.183),
                  g_B = c(-0.472,-0.589,-0.422,-0.479,-0.371,-0.112,-0.03,-0.134,-0.132,-0.12,-0.1,-0.227,-0.137,-0.025,-0.02,-0.011,-0.024,-0.377,-0.494,-0.49),
                  g_C = c(-0.278,-0.385,-0.31,-0.3,-0.279,-0.033,0.089,0.154,0.231,0.308,0.341,0.4,0.397, 0.383,0.37,0.309,0.248, 0.271,-0.212,-0.3),
                  g_D = c(-0.544,-0.604,-0.505,-0.541,-0.475,-0.173,-0.054,-0.026,0.098,0.199,0.257,0.338,0.37,0.372,0.36,0.274,0.217,-0.007,-0.552,-0.528),
                  g_E = c( 0.213,0.239,0.277,0.304,0.384,0.47,0.537,0.547,0.5,0.489,0.452,0.469,0.439,0.419,0.408,0.372,0.049,0.335,0.304,0.282),
                  g_F = c(0.234,0.248,0.3,0.324,0.391,0.444,0.512,0.503,0.448,0.434,0.407,0.421,0.4,0.385,0.376,0.335,0.087,0.323,0.295,0.275))

# Convert to the proper format
df2.melted <- melt(df2)

# Order of the factor in x-axis
df2.melted$Time <- factor(df2.melted$Time, levels = c("28Jan","13Feb","16Mar","1Apr","17Apr","3May","19May","4Jun","20Jun","6Jul","22Jul","7Aug","23Aug","8Sep","24Sep","10Oct","26Oct","11Nov","13Dec","29Dec"))

# Note group = variable to connect the points of the different groups
ggplot(df2.melted, aes(x = Time, y = value, colour = variable, group = variable)) +
  geom_line() +
  ylab("NDVI") + xlab("Time")

Very cool answer and succinct code

Thank you I really appreciate your help.

You are welcome @Sharbanoo! If the answer was useful to you please consider to accept it


text:
Empty values in FILTER field after running Mutect2
I was expecting something like in this image:

But when I ran Mutect2, the FILTER field was empty in the output VCF file. How can I enable PASS filters?
chr1    19063526    0   T   TAGCTTAGAATGGGATGTC 0   0   DP=1;ECNT=1;MBQ=0,34;MFRL=0,193;MMQ=60,52;MPOS=1;POPAF=7.30;TLOD=4.20   GT:AD:AF:DP:F1R2:F2R1:SB    0/1:0,1:0.667:1:0,0:0,1:0,0,0,1

Command used: gatk Mutect2 -R hg38.fa -I test.bam -O test.vcf.gz
GATK -v4.1.2.0
Mutect2 version 2.2

Linux v4.15.0-1044

Please [edit] your question and add i) the exact command you ran on the file and iii) the version of the tool(s) you used and the operating system you ran them on and iii) a variant line from the output VCF file.

@terdon It was edited, and now it includes what you asked for.

What is the **FILTER**? Please show us an actual line from your file. If the FILTER field was emtpy, then leave it empty. We need to have a representative line from your file that shows the problem.

chr1 19063526 0 T TAGCTTAGAATGGGATGTC 0 0 DP=1;ECNT=1;MBQ=0,34;MFRL=0,193;MMQ=60,52;MPOS=1;POPAF=7.30;TLOD=4.20 GT:AD:AF:DP:F1R2:F2R1:SB 0/1:0,1:0.667:1:0,0:0,1:0,0,0,1
is the actual line

Ah, OK thanks. So both FILTER and SCORE fields were 0, then? Not empty?

https://software.broadinstitute.org/gatk/documentation/tooldocs/4.1.2.0/org_broadinstitute_hellbender_tools_walkers_mutect_Mutect2.php      This link has tumor-only mode option.

Yes, sorry, my bad. I was looking at an older version. But please answer my previous question: were the FILTER fields empty or were they set to 0 as you show in your example?

I didn't set to 0. The output itself is 0. I tried --RF [PASS] even then i get 0 as output for all variants.


text:
Principle of TMT Tags in Multiplex Proteomics
I am new to proteomics research and analyzing mass-spec data. I am trying to learn more about the theory/principle of TMT tagging (a type of isobaric labeling) coupled with LC-MS/MS spectrometry. I understand that each tag has a
1) protein binding group (amine reactive ONLY?)
2) reporter group (so we know the tag came from which sample?)
3) normalisation group (to ensure that each tag has same mass?)
4) cleavage group (so tag can be cleaved from the peptide?)
Did I get the purpose of all the four parts in the tag?
Also, what does it mean by 4-plex, 10-plex, n-plex tag? And why do the tags have different isotopes?
If I could have some clarity on the workflow of TMT tagging and how it leads to analysis, it would be really helpful.

TMT stands for tandem-mass-tag. The idea is to tag different samples with one of the tags then pool all samples together and run them through the mass spectrometer together. Then the reporter groups are split off and quantified (ideally) for each peptide individually. There are different methods with different numbers of isobaric labels. Isobaric means all the labels have the same mass therefore they behave exactly the same in the mass spectrometer at first. E.g. TMT11 means you have 11 different tags with identical masses, so you can pool and process 11 samples at a time.
The protein binding group attaches the label to the peptide. Then the tagged peptides from all samples are run through the chromatography and the spectrometer. At first, the MS instruments simply collects the tagged peptides with different m/z ratios at each point in time for a while. When enough peptides are collected reporter groups are cut off and the peptides then quantified. This cutting-off happens by increasing the internal energy of the collected peptides. Because the cleavage group has the weakest chemical bond it brakes first and releases the reporter group. Again, this step separates the reporter group from the rest. Finally, the cut-off reporter groups are quantified.
The different reporter groups don't have the same mass, otherwise you would not be able to differentiate between them in the MS instrument. But you need same masses before quantification to ensure that the same peptides are all collected together. Therefore, you need to balance out the different masses of the reporter groups. That is the job of the normalization group. So, the normalization group ensures that overall the tags have the same mass and the peptides with different tags are not treated differently in the MS instrument before quantification. Then for quantification the reporter groups (with different masses between each other) are cut off and the MS instrument can therefore detect and quantify them individually. There is a time dependency of the whole process which happens in multiple phases. The whole process makes it possible to compare the amounts of peptides coming from the different samples.
The chemistry in the chromatography column is quite complicated which makes it difficult to compare samples from different experiments. The presence of one peptide might enhance or suppress the signal of another peptide or impact the chromatographic retention time. Therefore, you should only compare within the same sample. More tags TMT11, TMT16 etc. allow you to run more samples in one run.
You can then use one (or more) channels (mass-tags), to tag a standard sample and then compare everything to this standard that is assumed to be identical over multiple experiments. This way you can compare even more samples than available tags. e.g. more than 11.


text:
Error Viewing file on HiGlass
I am implementing hi-glass python on Jupyter, and upon running the test data code
from higlass.client import View, Track
from higlass.tilesets import cooler

ts1 = cooler('../data/Dixon2012-J1-NcoI-R1-filtered.100kb.multires.cool')
tr1 = Track('heatmap', tileset=ts1)
view1 = View([tr1])
display, server, viewconf = higlass.display([view1])

display

I get the error message stating:
Error Retrieving tileset info: http://localhost:60501/spi/v1

Please advise on the correcting. The code above has been copied from the HiGlass webpage.

Hi, welcome to the site. Did you "set up a temporary server" as the guide suggests ? It doesn't provide how to do it, but presumably it should be easy to do

I'm one of the authors of HiGlass. Sorry you're getting that error. I just tried it out in a fresh environment and it worked for me. Just one thing is missing in the docs which is that you need to do import higlass.
Here's what I did to set up my environment:
conda create --name test-bi python=3.6
source activate test-bi

The rest is from the docs:
pip install jupyter higlass-python
jupyter nbextension install --py --sys-prefix --symlink higlass
jupyter nbextension enable --py --sys-prefix higlass
jupyter notebook

It looks like either you or somebody else added a bug report here: https://github.com/higlass/higlass-python/issues/29
I'll respond and continue the conversation there.
Also, there's a repository full of examples from the HiGlass SciPy presentation here: https://github.com/higlass/scipy19


text:
How to efficiently get human gene names from NCBI based on a large list of SNPs
I found a good answer related to my question here:
How to get a list of genes corresponding to the list of SNPs (rs ids)?
But it says about small number of SNPs. I want to get gene names based on approximately 12000 to 19000 SNPs. NCBI says about 10 read per second with API key. Is there any efficient way to get gene names for thousands of SNPs? Is that database available to store locally? Any advice please?

Yes, you can download the NCBI data locally, see this page.

Modify SNPs and gene annotations for your genome assembly of interest. This example shows how one might do this for hg19, for instance:

Download SNPs and convert to BED:
$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz \
    | gunzip -c \
    | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - \
    > hg19.snp151.bed

Filter this BED file for entries of interest:
$ grep -wFf snps-of-interest.txt hg19.snp151.bed | cut -f1-6 > hg19.snp151.filtered.bed6

Download gene annotations and convert to BED:
$ wget -O - ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
    | gunzip -c \
    | awk '($3 == "gene")' \
    | gtf2bed \
    > gencode.v19.genes.bed

Note: In the future, newer versions of Gencode annotations may be available for your reference assembly. Always check their site to confirm what is most recent.
Map genes to SNPs:
$ bedmap --echo --echo-map-id hg19.snp151.filtered.bed6 gencode.v19.genes.bed > snps_with_associated_gene_names.bed

A list of genes that overlap a SNP of interest will be in the seventh column of the output.
If this must be done within Python, you could use something like subprocess.call(...), where ... includes commands wrapped in bash -c.
I don't recommend this because of the need to escape quote marks and deal with Python's quirks in running command-line tasks, which makes a Python-based approach very difficult to set up and maintain. I'd suggest learning some basic bash scripting to make this easy and fast.

I have SNP list. Can your suggested technique be applied in Python?

Maybe look into subprocess.call() to run command line tasks inside Python.

Run them through the Ensembl VEP. You'll get the gene name and all the information about how the genes are affected by the variants.


text:
How to predict protein binding to ligands?
I want to predict from protein sequence if the protein binds to metal, nuclear or small ligand. How can I do this ?  Which features are relevant if I want to use them in a machine learning algorithm ?

Welcome to the site! Could you please explain what features/data do you have available? Also it could be helpful too if you provided which machine learning algorithm do you want to try. Also could you clarify what does it mean "metal ligand" ? I've found this from wikipedia, but I doubt that this is what you mean

If you are training a machine learning algorithm you probably want to train against things you know to be true, rather than making predictions from the sequence and training against the predictions.
In this case you will want to use protein structures with the relevant ligands bound as your data. If you go to the RCSB PDB advanced search you can search for entries with any ligand present, specific ligands present, protein and DNA present, etc. Then you can extract the FASTA sequences for these entries.
See also databases like MetalPDB for specific ligand types, and one of our papers for an example of creating a dataset of metal binding proteins.
Features to use would be the raw sequence as a one-hot encoding on amino acid types. You could add other sequence level features such as secondary structure predictions, disorder predictions or physico-chemical properties. However a decent network would be able to predict much of this if it needs to. You could create alignments and put in a sequence profile or covariation data, which would be giving the network more information.


text:
Is the genome a programming language (i.e. LISP)? Can we analyze it with computer science theory?
Question moved here from the biology stack exchange.
It looks like the genome is a kind of LISP language.
Operon structure of genome:
- https://github.com/philschatz/microbiology-book/blob/master/resources/OSC_Microbio_11_07_Operon.jpg
LISP function structure:
- http://support.ircam.fr/docs/om/om6-manual/res/listprefix.png
Nested genes:
- https://player.slideplayer.com/31/9782148/data/images/img11.jpg
Nested LISP cons cells:
- https://www.tutorialspoint.com/lisp/images/treestructure.jpg
The basic similarity I am seeing is of general structure: (func var0 var1 ... varn) and nested cons cells.
The correspondence with genome structure is the operon portion, which has:

operator = function
structural genes = variables
structural genes can be nested like cons cells

Is this a global structure for the genome, then a specific prediction is we can do something simple like counting open and close markers in the nested genes and they will be equal.
So, why would such a correspondence matter, if true?  The basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming. Or just identify simple properties like balanced opening closing markers. The mere fact it has a LISP structure would be fascinating if true.  And since LISP is essentially an abstract syntax tree (AST), which all programming languages have, perhaps this structure I see in the genome is actually an AST, which implies an even more fundamental correspondence between the genome and programming languages.
More generally, we have a wealth of knowledge reverse engineering human programming languages.  If the genome is also a programming language similar in some respects to human programming languages, such as being a LISP variant, then we could apply this same reverse engineering knowledge base to better understand the genome.
So, first of all, is my high level observation correct?  Does the genome, at least for bacteria, have a global LISP like structure?  If yes, then if the genome is a kind of programming language, i.e. LISP, can we analyze it like software to understand its functionality?  Has anyone done this before?  Simple searching on Google Scholar turns up lots of software for analyzing genomes, but nothing about analyzing the genome as software.  It would be interesting to know whether anyone has seriously examined the question before, and either made some useful discoveries, or found out there is only a very loose analogy with programming languages, and thus knowledge regarding programming languages is not of use regarding the genome.
Hopefully this question is clear enough now, if not please ask me for further clarification.  I would like to get some sort of decent answer.
UPDATE: Someone who has put thought into "DNA is code" idea.  At a more general level than whether specific genome structures match specific programming structures.
UPDATE #2: The late Scott Federhen of NCBI wrote an article comparing genomic replication with the lambda operator.
UPDATE #3: It might be the case that LISP was inspired by the genome.

I'm afraid you're going to get much the same reception here as you did in Biology. Please [edit] and point to specific features of the genome that, to you, are similar to specific features of LISP. What would the "func" be in the genome? What are the vars? The genome is usually thought of in terms of natural languages with words (codons) making sentences (genes) and punctuation (start and stop codons). And I fear you are vastly underestimating the complexity of the cellular machinery that reads the genome. Also, your first link is broken.

AKA Richard Dawkins without the proofs. You need to draw down a specific testable hypothesis, the days when rhetoric ruled evolutionary theory ended with sequencing. LISP has a lot of fans (Perl people like it), biological systems are far less tractable and much more diverse.

https://xkcd.com/793/

@terdon I explained exactly what the correspondence is and a simple prediction that results.

Well, operons are only relevant to bacterial genomes. What do you mean by "structural" genes? And yes, genes can be nested, but they can also overlap on opposite strands, are spliced, up and down regulated, all sorts of things. And then there are exceptions, edge cases, outliers... I still don't see what you're getting at here to be honest.

@terdon it looks like the genome uses a general structure that is the same as the basic LISP structure, so any theory regarding LISP type languages, and programming languages in general, is applicable to understanding the genome.  I'm asking if my understanding is correct, and if anyone has undertaken this sort of analysis before, i.e. analyzing the genome as if it is a programming language.  Let me know if this is still unclear and what more I can do to clarify my question.

I guess it may come down to my own ignorance (I am a biologist, not a computer scientist) but I have no idea what kind of theory you would be able to apply or how. My only contact with LISP is my .emacsrc file. But many of us here aren't engineers, so your question would benefit from a more specific explanation of what sort of tools/approaches you are thinking of. More to the point, I still feel you are not grasping just how complex and, most importantly, stochastic and chaotic system a genome is. I just don't see what benefit your approach would bring or how it would even be possible.

@terdon my basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming.  Or just identify simple properties like balanced opening closing markers.  The mere fact it has a LISP structure would be fascinating if true.

@terdon also, more generally, I'm asking if anyone has taken this approach before, and demonstrated its effectiveness or found out there is some major problem so analyzing the genome as a programming language does not yield anything useful for a particular reason.  Do you happen to know either way what research has been done in this area?  Intuitively it seems a useful idea, because there is a wealth of knowledge regarding reverse engineering of programming languages, so we might be able to apply this knowledge to reverse engineer the genome.

Please [edit] your question and explain all this. Ideally, explain why you think programming languages are a good paradigm and better than natural languages. I still haven't seen any good evidence presented here and fear it might be a case of "when all you have is a hammer, everything looks like a nail". Granted, perhaps I'm the one with the hammer and not you, but that's why it would be really helpful if you could express your idea in terms of specific features of programming languages and how those would help you understand the genome better than we do already.

@terdon I have added the explanation from the comments to the question.  Hopefully it is clear enough now.

I think the name we use for "reverse engineering knowledge base to better understand the genome" is genome annotation.
But more seriously. I don't think there is any person who would possess sufficient knowledge of both to find practical analogies. Because we have plenty of ways how to analyze genome from different perspectives (e.g. genomes as subjects of mutations and selection over millions of years)m hence trying to fit some cs theories might be... difficult.

@KamilSJaron surely one more option might be helpful?  Why do you think this question is encountering so much resistance?  I've tried clarifying it a number of times, and in my opinion it is pretty clear what I am asking.

I'm voting to close this question as off-topic because it has nothing to do with bioinformatics

@MichaelG. seems like interpreting the genome as code is precisely bioinformatics.  Per the Wikipedia definition: As an interdisciplinary field of science, bioinformatics combines biology, computer science, information engineering, mathematics and statistics to analyze and interpret biological data.

@yters: bioinformatics is interpreting biological data with computers. The opposite is sometimes called "computational biology" which is using biological systems AS computers. Or analysing biological systems as if they were a computer. So DNA computing, for example.

@gilleain this is a bit different than using the systems as computers.  It is interpreting the systems as IT constructs, in this case a programming language.  Reading the Wikipedia articles on both topics, it still seems to me that bioinformatics is the best match, since the idea I propose relies on using computer science concepts to analyze the genome.  There is at least a surface visual similarity that seems quite strong, so I'm curious whether people have looked into this idea and what they have discovered.

I see what the metaphor has been inspired by, where I think it got turned around, and might have a better one.
I'd like to back up a bit and establish some common ground.  It isn't wrong to expect that some aspects of genetic code and software code are similar, and while biology is usually very specific, some aspects can be seen productively from the light of computer science ideas.  For example, mass-action molecular dynamics, which are often a good model for protein interactions, are Turing-complete.  Here's a fun example where people have written a programming language out of molecular dynamics:  https://arxiv.org/abs/1809.07430
It's also true that proteins have particular substructures leading to certain functions, such as using alpha-helices to react less with the outside world.  Learning about larger substructures, called protein domains, help identify the functions of families of proteins.
The problem is that the parts of the genome which appear nested like LISP cells actually aren't included in the protein.  The regulatory regions, promoters, inhibitors, introns, and so forth never become part of the key "computational" machinery.  This nesting isn't actively part of the computational activity proteins do perform.  The main quality of LISP is that the structure and function are perfectly self-similar throughout an expression, but the structure of the regulatory regions are neither replicated down into the protein nor are necessarily recapitulated on a larger scale.
Unfortunately, we don't expect the computational structure of proteins to be the same from one protein to another.  This is because "binds only with X" and "binds only with Y" must be implemented differently to successfully be distinct operations, and it is these binding constraints by which molecular dynamics implement "computational" operations, such as "not X" and "not Y", and so those are distinct as well.
I think a better way to think about the regular structures of the regulatory apparatus is more like a network protocol.  Along an otherwise undifferentiated binary channel, one has to identify when and where the data is delivered to an application.  Similarly, along a strand of bases, the machinery of translation has to have an indication of where the "application" of the protein starts and ends, and when that protein should be delivered to the cell.  We might compare some of the images of genetic structure that provided in the original question with layouts of network headers: http://slideplayer.com/slide/4798836/15/images/16/TCP/IP+Packet+Structures.jpg  We might compare bacterial versus multi-cellular regulation as having different delivery mechanisms and data structures for delivering the application information, for different functional consequences, in the same way TCP and UDP have different structures for different trade-offs.
This is all complicated somewhat by the fact that proteins do then interact with the regulation of genes, leading to interesting "computational" dynamics such as feed-forward regulation.  However, I'm not aware of any particular structure of the regulatory apparatus that necessarily correlates with these larger functional patterns.  The opposite is frequently biologically useful: the same structure of a protein helping to activate one gene might simultaneously inhibit another.
One book I've found very accessible and useful in understanding where computation is useful biologically is Dennis Bray's "Wetware: a computer in every living cell".  It gave me a clearer picture of how computation occurs in cells.  That should open up some introductions to systems biology which explains the biological-role of larger molecular "computational" patterns.

Great answer!  First answer of substance on either site!  I'll have to think it over, and if I find your reasoning persuasive, I'll accept this answer.

Thanks!  I'm glad you found the answer constructive.  After I found a useful perspective, I enjoyed thinking about the the question.

You say one problem with the idea is that the LISP like structure doesn't actually become part of the computational machinery.  I'm not sure this invalidates the connection I see.  When a LISP program is compiled down to the machine code format, the high level structure disappears.  This sounds like what you are saying with regard to proteins.  All the structural patterns disappear once the genome is compiled into a protein.  So, to my ears it seems you are saying LISP and the genome are the same in this regard, instead of different.  Perhaps you can clarify?  BTW, good point about protocols.

A high level comment, I think your network protocol idea and this LISP idea are consistent with a broader idea that the genome may have a number of different protocols that are followed, just like a modern distributed computing system.  And while they are understood holistically in reference to each other, each protocol can also be analyzed by itself.  Perhaps the same mode of thinking applies to the genetic code and regulation networks.  While they can be understood holistically, they perhaps can also be understood individually.

Right, so it is true that a kind of compilation does happen, where sequences of bases undergo a translation into the amino acids that make the protein.  The observation that I'm making is that the part of the gene in which you observed the LISP structure, is not included in that translation process.  Instead, those structures have been stripped away at that point.  This is similar to how when data is provided by a network protocol to the application, those headers are not delivered to the application.

The question is then whether or not we expect to see regular structures corresponding to computational structures in the part of the gene that is then translated, such as how we would see machine instructions for branching, adding, and these sorts of things.  Unfortunately, the answer is generally not, as the computations implied have to be specific for each different input.  If the answer were yes, the challenges of bioinformatics would be much easier!

It occurs to me that it might be helpful to know that promoters, inhibitors, introns, and so forth are occasionally called "non-coding regions" given that they aren't included in translation.  The following diagram might be helpful:  https://en.wikipedia.org/wiki/File:Gene_structure_prokaryote_2_annotated.svg

I see, just the genes inside the segment remain when translation takes place.  So, in that sense it is dissimilar to LISP, since the function in a LISP s-exp would be a part of the compilation process.  On the other hand, the non-coding section does seem to play a role in whether and how the genes are translated, if I'm understanding the answers so far.  And, thus this prokaryote genome structure still looks a lot like a LISP s-exp, since it can presumably also be nested: https://en.wikipedia.org/wiki/S-expression

As you've inferred, it is entirely right that these non-coding sections do play a role in regulation.  The way that regulation generally works, crudely speaking, is that whether transcription can go ahead or not depends upon how proteins bind to these regulatory sections, either guiding or blocking the transcription machinery.  The content of these regulatory sections must then be specific to the proteins that do work, so that once again the computational work that is done unfortunately must be specific to each family of proteins and not part of a more general computational mechanism.

I agree that there is a kind of nesting that can be represented as an S-expression.  Instead, it's my contention that the S-expression could be a misleading choice for representing that structure, as I propose that the nesting that you've observed is better explained as the nesting of protocols rather than the recursive nesting of context-free grammars.  By protocol nesting, I mean that after the data has been delivered from one protocol, and the metadata of that protocol is removed, the resulting data has the data and metadata of the next protocol.

Thanks again, this is a truly excellent series of answers, and you've clarified the issue greatly in my mind and given me a reference to further my understanding.  Great work!

Great, I'm glad to read that this has been helpful.  I have to hand it to you for not being dissuaded from pursuing a potential hypothesis for the structure you observed, while yet entertaining suggestions counter to your hypothesis, until you found an explanation that attempted to address all the issues.

So, first of all, is my high level observation correct? Does the genome, at least for bacteria, have a global LISP like structure?

No. I'm afraid this is what it boils down to. I don't see any real similarity to LISP or other programming languages at all. Perhaps it would be clearer if you could explain in detail why you think they're similar, but as it stands, no, there is no such structure apart from a very simplistic similarity which is sort of the same as saying that the genome is like a car because it carries genes and could therefore be studied by a car mechanic.

You don't find the correspondence list in my question clear enough?  I'm asking about the structure, which seems to be the same as a LISP program, or perhaps more generally an abstract syntax tree.  And yes, the structure itself is pretty simple, but still there seems to be an exact correspondence, more so than referring to a car. I am somewhat bemused this question is so hard for people to understand and answer.

@yters no,  I don't find it clear at all. As far as I can tell, you seem to be comparing functions and their arguments to operons and their genes, but I don't see any similarity at all there. That's why I posted so many comments asking you to get more specific. I have no idea why you feel that operons are more like functions than they are like cars or bowls or anything else that can contain items. You might want to look into gene regulatory networks though, those are closer to what you're getting at.

You could say there is a disanalogy in that the promoter/operator is not really a function, but controls whether the genes in the operon are "executed".  But there is still generally a notion of "execute this thing" and a nested tree structure.

@yters not really. That's not how the promoter works. Promoters tend to control how much it is executed, not whether. And they're just the tip of the iceberg in terms of regulation.

but is it correct the promoter/operator preamble somehow controls the manner in which the subsequent genes are 'executed'?

@yters it is only one of many things that affect gene expression, some of which are before the gene, some in it, some after it, some can even be quite distant from it and others aren't even part of the genome at all. As I said before, you are vastly underestimating the complexity of the genome.

I don't deny the extreme complexity of the greater context.  But, just as I can understand the C code of some UNIX utility without having to understand UNIX, perhaps we can analyze a genome by itself, to at least some degree.

You mention this as the base similarity you find:

operator = function
       structural genes = variables
       structural genes can be nested like cons cells

As other has pointed, this just considers the genome as a sequence of where only the linear order is important. It doesn't take into account other things that happen on the genome that are know to affect the cells like the chromatin structure, the duplication of certain regions, the transversion, the effect of moving elements, methylation, or simply that a single gene can have multiple functions (pleiotropism, moonlighting) while a programming function does not.
Even if we considered it as code, certainly by using tools like HMM we could also learn more about the structure of the genes we have. While certainly, we need more tools than just the sequence to discover what does each gene do.
There are a lot more differences than similarities on this comparison, it can be helpful for certain aspects but it fails for other aspects of cells (without even going into the virus discussion or the mitochondrial).

The only possible parallel was the "processor 'core' wars" in Avida and Tierra. Tierra was initiated by Tom Ray, around 1990. Here they had a 'population' of compiled binaries were randomly mutated (yep the binary code was mutated) and reproduction was dependent on how much of the processor could be consumed. The simulation worked and an in silico ecosystem evolved, but was heavily dependent on the preconditions of the simulation and proved no more reliabe at understanding evolution than differential equestions. So it was a nice idea, and very cool it worked, but it came to zero.
That aside the best example was the concept of the "selfish gene", but has very little to do with programming languages IMO and everything to do with understanding kinselection.

Not sure what question you are answering here.  I am not asking about evolution simulations.

To say I want "genome understanding" without evolution, is like saying I want a watch without hands. Can't have one without the other.

I would say it is more like I can understand the C source code of a UNIX utility without having to know the history of UNIX or C's development.

Believe it or not, there are people who know more than you do about programming, and who also know a hell of a lot more than you do about the genome.  This is not an avenue those people think is worth pursuing.
DNA is a chemical, and interacts with other chemicals based on its shape.  There is no equivalent to that in a human designed programming language.

This seems unnecessarily combative. The OP never claimed to be an expert!

Which is why OPs first step should have been to read what experts in both subjects have been doing in the field he's interested in, instead of trying to impress a bunch of volunteers with programming jargon.

Wow. Please take a step back and stop assuming malice. The OP did do some research and did try and see if any work had been done. Nobody's out to impress anybody here, so please stop being so aggressive!


text:
swap genotypes for specific fields of a vcf file
Is there a tool already established to "swap" the ref/alt information in the GT, AD and PL fields of a vcf file?
For the GT field, I need to transform: 0/0 -> 1/1 ;  1/1 -> 0/0 ;  0/1 -> 1/0 ; and finally 1/0 -> 0/1
For the AD field, I need to swap the columns: x,y -> y,x
For the PL field, I need to swap the columns: x,y,z -> z,y,x
Sample input:
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  0/0:17,0:17:48:0,48,720 0/1:14,8:22:99:167,0,347 .....

Sample output:
    Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:720,48,0 1/0:8,14:22:99:347,0,167 .....

note: the same operations will need to be carried out for all samples if the VCF file has multiple samples (the ... stand for other potential samples).
Before I dive into scripting ( basically playing with regex), I would like to know if something to process this is already available. I looked into bcftools and fixref but they don't seem to help for this.

What do you mean by "note: the ... stand for other columns to be swapped as well,"? You mean carry out the same operations on the GT, AD and PL of any other samples in the VCF file?

yes, better wording here! thanks

How generic does this need to be? I mean, is it enough to give you a solution that works specifically for the order of fields you have (GT:AD:DP:GQ:PL) or do we also need to be able to cope with different orders (so need to infer the position of the AD field from the INFO field)?

for now, a solution for the order of fields I have is okay. I don't have a lot of experience with vcf file, and I thought that the order of fields was "fixed" for all vcfs :9

If we can always assume your fields will be GT:AD:DP:GQ:PL, you can use this relatively simple, brute approach awk script:
awk -vOFS="\t" -F"\t" '{
  ## just print the header lines
  if(/^#/){
    print;
  }
  else{
    ## The number of samples will be the extra number of fields
    ## in addition to the default VCF 9
    sampleNum = NF - 9
    ## Fix each sample
    for(i=1; i<=sampleNum; i++){
      ## get the fields
      num=split($(9+i), fields,":");
      ## genotype
      if(fields[1] ~ "0/0"){ fields[1] = "1/1"}
      else if(fields[1] == "1/1"){ fields[1] = "0/0" }
      else if(fields[1] == "0/1"){ fields[1] = "1/0" }
      else if(fields[1] == "1/0"){ fields[1] = "0/1" }
      ## AD
      split(fields[2],values,",")
      fields[2] = values[2]","values[1]
      ## PL
      split(fields[2],values,",")
      ## Now set these as the sample data
      $(9+i) = fields[1];
      for(k=2;k<=num;k++){
        $(9+i) = $(9+i)":"fields[k]
      }
    }
    print;
  }
}' file.vcf

On your example input, this returns:
Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:0,48,720 1/0:8,14:22:99:167,0,347    .....

thanks @terdon, I guess this means we don't really have a software able to handle it and I'll need to dive into scripting :) your solution works so i'll accept it, although I'd prefer to use some available software in this case.

@aechchiki yeah, I doubt there will be software for this since you are doing something that is very much not standard: you are swapping data around between different alleles.


text:
how to remove range from fasta header
Could you please suggest me how I can remove range from fasta header:like these number from below sequences which has some colons indicating the range of the genome
:147010-147657
:149201-149845
> FVJT01000004.1:147010-147657 Mycobacterium abscessus subsp. massiliense strain 974
GTGACCAGCAGCGCAGACCGATCCCCCGTCGAGGCGATGGTGAGCCACGCCGAGGGCGCCATCTCCGAGG
ATGACATCGTCGCCGCCGCCCGGGAACGTGCCGTTGACCTGGGCGCGGGCCCGGTGACACCCGCGGTAGG

>FSAX01000002.1:149201-149845 Mycobacterium abscessus subsp. bolletii strain 286 genome
ATGACGACGCCCGTCAGCTTCGATAACAGATTCTTCTCCGCAACGTGGTCTGCTATCGCGCGGCGCGAAC
CTCAGGAGATCAAGGAGCTGCGTGCGGAGAATCTACGCGGCCTGACCGGCCACGTTCTGGAGGTGGGAGC

Invest some time in learning unix tools such as grep, sed and awk.
sed -r '/^>/s/:[0-9]+-[0-9]+//' input.fasta > output.fasta

This line of code matches the pattern :<numbers>-<numbers> on lines that start with > (FASTA sequence headers) and replaces that pattern with '', essentially removing text that matches the pattern.

Thanks Ram, could you please suggest any book related to unix command line that can explain all these small problems related to bioinformatics.

I'd recommend looking at Rosalind, but I think that's mostly python. Unix utilities are learnt as you stumble upon problems, as there's seldom some sort of comprehensive exercise to learn these utilities. Try googling "learn sed" or "learn grep" or "learn awk" and following tutorials that pop up. For awk, you can start with an Excel sheet of data, export the data to a tab-separated text file and use awk to do things you'd normally use Excel filters to do.

Alternatively, Perl provides a slightly simpler regular expression syntax (\d instead of [0-9]) that I prefer over sed.
perl -pe 's/:\d+-\d+//' input.fasta > output.fasta

Using -pe is even simpler: perl -pe 's/^(>\S+):\d+-\d+/$1/' input.fasta > output.fasta


text:
Produce a single sequential FASTA sequence out of BAM
I'm having problems properly looking for a solution because I'm a layman in Bioinformatics not familiar with the terminology. I'm hoping you can nudge me in the right direction, please! Thank you very much.
What I want in the end
A FASTA file containing a single sequential sequence (chromosome 1-22..) of the human genome.
What I have

23gb of "clean_data" which are .fq.gz files.
21gb of "result_alignment" which is .bam+.bai
"result_variation" indel.vcf.gz and snp.vcf.gz

What I tried
I used samtools fasta reads.bam > reads.fasta which results in many lines like so:
>V100005418L4C001R001322578/1
AAATAGCCCACAC...
>V100005418L4C001R006084507/1
GCCATAAAGCCTA...
...

But these include many many duplicate lines and the reference genome from the human genome project GRCh38 is only 3GB in size, so I thought that this FASTA sequence needs to be properly "aligned" somehow?
Am I thinking correctly here?
I'm really confused with all of the terminology and don't know exactly what I need to look for. Could someone kindly put the process into simpler words for me?

What you have in BAM format is an alignment of reads to a reference. What you are looking for (a single fasta per chromosome) is a new assembly. Using "samtools fasta" will just get you each read in fasta format, which is clearly not what you want. In addition to doing a (de novo) assembly of your reads you could make a (reference-guided) consensus sequence. You can find some pointers on how to do that here: https://www.biostars.org/p/367960/

Thank you so much for your comment @Wouter. The link you posted sounds very promising. I will try it out and get back to you. Meanwhile I clarified in my post what kind of data I have. Maybe you want to take a deeper look. Thanks again!

another option is to take the reference genome you used for alignment, and 'mutate' it with the variants you identified. One tool to do that is GATK FastaAlternateReferenceMaker. I don't know how well it deals with indels though.
And keep in mind that the sample which you sequenced is diploid. If you want one fasta sequence per chromosome for that individual then you have to think about how to handle heterozygous sites. If you want one paternal and one maternal fasta then you are in bigger trouble because you need to know on which phase the variants are relative to each other.

What I wanted to do is called "Consensus Sequence". Two steps were needed:
$ bcftools mpileup -Ou -f ref.fa input.bam | bcftools call -Ou -mv | bcftools norm -f ref.fa Oz -o output.vcf.gz

After that you can create the consensus:
$ bcftools consensus -f ref.fa output.vcf.gz > out.fa

If you just want the human chromosomes in a FASTA format, why not downloading directly the individual chromosome files (chr*.fa.gz)?
If you don't rely on the "official" releases (for whatever reason), then what you need is first to assemble your reads into genomes, possibly (but not surely) going down to a chromosome-level resolution. Different approaches for assembly are possible/available, mainly depending on (1) what do you biologically need to answer and (2) the type of sequencing data you have - we'll need more details to be able to help you in this direction.
The samtools fasta utility will "just" convert your reads (likely given from a sequencing facility) from an alignment (BAM) format into a reads (FASTA) format, there is no concept of chromosomes here. You won't get assembled chromosomes this way, just extract the sequence in FASTA format from your alignment file.

Thank you very much for your input! I don't have something biologically to answer, this project revolves about sequencing(?) this specific individuals DNA. I need the DNA from "start to finish" so to speak, Chr1 to Chr22, X,Y in ascending order, is what I'm trying to say. I understand now that the FASTA conversion is not what I need, I need to assemble my reads first. As for data I have: 23gb of "clean_data" which are .fq.gz files. 21gb of "result_alignment" which is .bam+.bai and I got "result_variation" indel.vcf.gz and snp.vcf.gz. Can I use these to achieve what I want? Thank you very much!

Usually when you sequence a person's genome, you're trying to identify variants that the individual has, compared to the reference genome. In that case, you should do what @aechchiki said, and download the official reference FASTAs, and create a VCF file containing all the variants for that individual.
Alternatively, it also sounds like you might be trying to do de novo genome assembly. That is a very different topic from variant calling, and that would give you a "reference" that is very different from the reference genome builds mentioned above

@DanielBiegler please [edit] your question and tell us why you are sequencing this individual since the data you will require will depend on what you are trying to do. It is extremely unlikely that you need the individual's full genome rather than the specific locations where their genome differs from the reference (variants). But if you don't us what biological question you are trying to answer using the individual's sequence, then we can't really help much.


text:
How to create Seurat object while RNA expression and ADT combined into one matrix
I got an input matrix which RNA expression and ADT capture are combined into one file. I loaded the file into Seurat successfully, however, when I tried to create Seurat object, it threw out an error saying
Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) :
  No cell names (colnames) names present in the input matrix

How can I create Seurat object when RNA expression and ADT data are combine together?

It shouldn't be a problem to have both types of data into a matrix as long as they are identified (i.e different naming of the columns). From the error that you say it seems that there are no column names in your matrix? Did you try to add colnames and see if it works?

@plat I checked the data and the columns(cell)  name is there. I don't know why it keep saying no cell names presents.

Then I don't know what is going on. Could you please edit your question with a minimal reproducible example (https://stackoverflow.com/help/minimal-reproducible-example) of your problem?

Given the fact that you have RNA and ADT data (probably on the same set of cells), some colnames (which represent individual cells) would be duplicated. Normally duplicated colnames are tolerated with matrices and sparse matrices, however, Seurat apparently does not do so, below is from the ?CreateAssayObject:
Details
Non-unique cell or feature names are not allowed. Please make unique before calling this function.

Possible solution to your problem: Seuart has a dedicated vignette for working with multimodal data and as you would see you will need to initiate your Seurat object with one matrix per assay: RNA and ADT. All you need to do is split your matrix into RNA and ADT, create your Seurat object with RNA data and then add the ADT data with:
seurat_obj_with_rna_only[["ADT"]] <- CreateAssayObject(counts = your_adt_matrix)
For efficiency, Seurat uses sparse matrices so don't forget to convert your data matrices to sparse.

I don't know if you solved the problem. Howerver,I think it's caused by multiple data types.
For output from CellRanger >= 3.0 with multiple data types
data <- Read10X(data.dir = data_dir)

This step creates multiple matrices, depending on how many types of data there are.
seurat_object = CreateSeuratObject(counts = data$`Gene Expression`)

So (counts = data$Gene Expression),You should also be able to choose the type of research you want to do.
I hope my answer will help some people.