text:
Shell script to validate fastq issue
I have two GSM ID as test case where both of them are having data when I checked through sra explorer tool but when I try to fetch through a script only one of them returns output where the other one fails.
GSM ID
GSM3603268
GSM2683458
Output
For GSM3603268
$ cat GSM3603268.test.fastq
@GSM3603268.1 1 length=100
TTTTTTAAAGAAAACTTGAGCTTTTGGAGCATGGCAACCTAGCCTGCAGACACCGTATCCCCTTGTCCACTTCCCCCTGCAAACCATAAGTCCATTCCTA
+GSM3603268.1 1 length=100
=@@DDDDA++2CDABCEG@GEHHEHGA9EDFHGHD3::B??<DBFFBFB===)5;CA;(7=?;?CE;7)7;>C((583;<9>(,89@@@@CCCCDD:BC:
@GSM3603268.1 1 length=100
AGTTCCACGAGTTTTCTTTTTTTTAAGTGGTAGGAATGGACTTATGGGTTGCAGGGGGAAGTGGACAAGGGGCTACGGTGTCTGCAGGCTTGGTTGCCAT
+GSM3603268.1 1 length=100
:?=BD:D?C)<CDGH4+<CCCGHI)??D<D?00?D*?39C<CG4BFE28;4(.7AB############################################

$ cat GSM3603268.test.log
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:18:46 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
Read 1 spots for GSM3603268
Written 1 spots for GSM3603268

Output for GSM2683458
$ cat GSM2683458.test.fastq
    bash-5.1#

$ cat GSM2683458.test.log
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:16:55 fastq-dump.3.0.1 int: string unexpected while executing query within virtual file system module - multiple response SRR URLs for the same service 's3'
2023-09-26T21:17:11 fastq-dump.3.0.1 warn: directory not found while opening manager within virtual file system module - can't open NC_000076.5 as a RefSeq or as a WGS
2023-09-26T21:17:11 fastq-dump.3.0.1 err: directory not found while opening manager within virtual file system module - failed GSM2683458

=============================================================
An error occurred during processing.
A report was generated into the file '/root/ncbi_error_report.txt'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.
=============================================================

fastq-dump quit with error code 3

The script
#!/bin/sh

set -x
PS4='[\\d \\t] '

# Check parameter for error
check=0
# Print fastq-dump executable path
echo $(which fastq-dump)

# Function to download FastQ with retries
download_with_retry() {
    local id="$1"
local retries=5
for attempt in $(seq "$retries"); do
    # Download start of fastq
    fastq-dump $(get_ngc) -X 1 -Z --split-spot "$id" > "${id}.test.fastq" 2> "${id}.test.log"
    numLines=$(cat "${id}.test.fastq" | wc -l)
    if [ $numLines -gt 0 ]; then
            echo "${id} has data... OK"
        return 0
    else
        echo "${id} does not have data on attempt $attempt... Retrying in 5 seconds..."
        sleep 5
    fi
done
echo "${id} could not be downloaded after $retries attempts... ERROR"
    check=1
    return 1
}

# Loop through all parameters to check validity
for file in "$@"; do
cp "${file}" .
    # Extract filename for sampleID
    file_basename=$(basename "${file}")
    id="${file_basename%".id"}"
# Start validation with retries
echo "Checking ${id}..."
    download_with_retry "$id"
done

# Exit with error if some FastQs are not accessible
if [ $check -gt 0 ]; then
    echo "ERROR: One or more samples have inaccessible FastQs.. exiting"
    exit 1
fi

What exactly is going wrong in the case of the other ID? No matter how many retries I give, it fails.
Any suggestion or help would be really appreciated.
The error log:
$ cat ncbi_error_report.txt
<Report>
 <Run>
  <Date>
   <Start value="Tue Sep 26 2023 9:16:55 PM"/>
   <End value="Tue Sep 26 2023 9:17:11 PM"/>
  </Date>
  <Home name="HOME" value="/root"/>
  <Cwd>/tmp</Cwd>
  <CommandLine argc="6">
   <Arg index="0" value="fastq-dump"/>
   <Arg index="1" value="-X"/>
   <Arg index="2" value="1"/>
   <Arg index="3" value="-Z"/>
   <Arg index="4" value="--split-spot"/>
   <Arg index="5" value="GSM2683458"/>
  </CommandLine>
  <Result rc="RC(rcVFS,rcMgr,rcOpening,rcDirectory,rcNotFound)"/>
  <User admin="true"/>
 </Run>
 <Configuration>
  <Files count="2">
   <File name="/etc/ncbi/settings.kfg"/>
   <File name="/root/.ncbi/user-settings.mkfg"/>
  </Files>
  <refseq state="not found"/>
  <krypto state="pwfile: not found"/>
  <sra>
   <quality_type>raw_scores</quality_type>
  </sra>
  <Config>
  <ConfigurationFiles>
/etc/ncbi/settings.kfg
/root/.ncbi/user-settings.mkfg
    </ConfigurationFiles>
    <APPNAME>"fastq-dump"</APPNAME>
    <APPPATH>"/tmp/"</APPPATH>
    <BUILD>"RELEASE"</BUILD>
    <HOME>"/root"</HOME>
    <HOST></HOST>
    <LIBS>
      <GUID>"119c217a-7b81-47e8-91d6-56d19c8c9f15"</GUID>
      <IMAGE_GUID>"119c217a-7b81-47e8-91d6-62229c64ee59"</IMAGE_GUID>
    </LIBS>
    <NCBI_HOME>"/root/.ncbi"</NCBI_HOME>
    <NCBI_SETTINGS>"/root/.ncbi/user-settings.mkfg"</NCBI_SETTINGS>
    <OS>"linux"</OS>
    <PWD>"/tmp"</PWD>
    <USER></USER>
    <VDB_CONFIG></VDB_CONFIG>
    <VDB_ROOT></VDB_ROOT>
    <kfg>
      <arch>
        <bits>"64"</bits>
        <name>"56d19c8c9f15"</name>
      </arch>
      <dir>"/root/.ncbi"</dir>
      <name>"user-settings.mkfg"</name>
    </kfg>
    <libs>
      <cloud>
        <report_instance_identity>"true"</report_instance_identity>
      </cloud>
    </libs>
    <repository>
      <user>
        <ad>
          <public>
            <apps>
              <file>
                <volumes>
                  <flat></flat>
                  <flatAd>"."</flatAd>
                </volumes>
              </file>
              <refseq>
                <volumes>
                  <refseqAd>"."</refseqAd>
                </volumes>
              </refseq>
              <sra>
                <volumes>
                  <sraAd>"."</sraAd>
                </volumes>
              </sra>
              <sraPileup>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraPileup>
              <sraRealign>
                <volumes>
                  <ad>"."</ad>
                </volumes>
              </sraRealign>
              <wgs>
                <volumes>
                  <wgsAd>"."</wgsAd>
                </volumes>
              </wgs>
            </apps>
            <root>"."</root>
          </public>
        </ad>
      </user>
    </repository>
    <sra>
      <quality_type>"raw_scores"</quality_type>
    </sra>
    <vdb>
      <lib>
        <paths>
          <kfg>"/usr/local/bin"</kfg>
        </paths>
      </lib>
    </vdb>
  </Config>
  <RemoteAccess available="false"/>
  <CurrentProtectedRepository found="false"/>
 </Configuration>
 <Object path="https://sra-pub-run-odp.s3.amazonaws.com/sra/SRR5755657/SRR5755657" type="database" fs_type="unexpected">
  <Dependencies>
   <List count="22" missing="22">
    <Dependency index="0" seq_id="NC_000067.5" local="false" path=""/>
    <Dependency index="1" seq_id="NC_000068.6" local="false" path=""/>
    <Dependency index="2" seq_id="NC_000069.5" local="false" path=""/>
    <Dependency index="3" seq_id="NC_000070.5" local="false" path=""/>
    <Dependency index="4" seq_id="NC_000071.5" local="false" path=""/>
    <Dependency index="5" seq_id="NC_000072.5" local="false" path=""/>
    <Dependency index="6" seq_id="NC_000073.5" local="false" path=""/>
    <Dependency index="7" seq_id="NC_000074.5" local="false" path=""/>
    <Dependency index="8" seq_id="NC_000075.5" local="false" path=""/>
    <Dependency index="9" seq_id="NC_000076.5" local="false" path=""/>
    <Dependency index="10" seq_id="NC_000077.5" local="false" path=""/>
    <Dependency index="11" seq_id="NC_000078.5" local="false" path=""/>
    <Dependency index="12" seq_id="NC_000079.5" local="false" path=""/>
    <Dependency index="13" seq_id="NC_000080.5" local="false" path=""/>
    <Dependency index="14" seq_id="NC_000081.5" local="false" path=""/>
    <Dependency index="15" seq_id="NC_000082.5" local="false" path=""/>
    <Dependency index="16" seq_id="NC_000083.5" local="false" path=""/>
    <Dependency index="17" seq_id="NC_000084.5" local="false" path=""/>
    <Dependency index="18" seq_id="NC_000085.5" local="false" path=""/>
    <Dependency index="19" seq_id="NC_000086.6" local="false" path=""/>
    <Dependency index="20" seq_id="NC_000087.6" local="false" path=""/>
    <Dependency index="21" seq_id="NC_005089.1" local="false" path=""/>
   </List>
  </Dependencies>
 </Object>
 <SOFTWARE>
  <VDBLibrary vers="2.7.47"/>
  <Build static="true">
   <Module name=""/>
  </Build>
  <Tool date="Nov 18 2022" name="fastq-dump" vers="3.0.1">
   <Binary path="/usr/local/bin/fastq-dump" type="alias" md5="c461c39bfa514aff3c4f7c0416ced617">
    <Alias resolved="fastq-dump.3">
     <Alias resolved="fastq-dump.3.0.1">
      <Alias resolved="sratools.3.0.1"/>
     </Alias>
    </Alias>
   </Binary>
  </Tool>
 </SOFTWARE>
 <Env>
 </Env>
</Report>

How are you running the script? What arguments are you giving it? What is get_ngc? What are the files you are copying (cp "${file}" .), what do they contain? And why are you running this as root!? Never run things as root unless it is absolutely necessary. Any minor error or bug could result in destroying your system if you run the wrong command as root!

to test it locally I i first start the container and then I copy the script inside the container and run using a single GSM id as input here the issue is for this particular series of GSM id none of them works

'and why are you running this as root!? Never run things as root unless it is absolutely necessary' im running it in a container

Use the SRR accessions, not the sample-IDs you are using.
With entrez e-utils, I found the SRR of GSM2683458:
esearch -db sra -query GSM2683458 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID is SRR5755657.
For GSM3603268, using esearch -db sra -query GSM3603268 | efetch -format xml | xtract -pattern IDENTIFIERS -if PRIMARY_ID -contains SRR -element PRIMARY_ID, the SRR is: SRR8571942
Then just do fasterq-dump SRR5755657; fasterq-dump SRR8571942.

well the pipeline is already made to use GSE id so i have no option but to use it


text:
How do I rewrite a read group using pysam?
I am trying to rewrite a SAM/BAM file with altered read group entries using pysam. In this simplified version, I want to take a BAM, and rewrite the SM tags in all read groups to the same string and copy the alignments with new header into a new BAM. The contents of the read group dictionaries I am changing appear to have no effect on the header that prints or is written to file output. How can I alter the read groups?
import argparse
import pysam

def rewrite_read_groups(input_filename, output_filename, sample):
    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        # replace header fields
        output_header = bam_input.header
        read_groups = output_header['RG']
        print(read_groups)
        for read_group in read_groups:
            read_group['SM'] = sample
            print(read_group)
        output_header.to_dict()
        print(output_header)
    finally:
        bam_input.close()

    print(output_header)

    try:
        bam_input = pysam.AlignmentFile(input_filename, 'rb')
        bam_output = pysam.AlignmentFile(output_filename, 'wb', header=output_header)

        # copy alignments
        for read in bam_input.fetch(until_eof=True):
            bam_output.write(read)

    finally:
        bam_input.close()
        bam_output.close()

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Rewrite bam read groups")
    parser.add_argument('-s', "--sample", help="Set sample for all read groups")

    parser.add_argument('-i', "--input_file", help="Input bam", required=True)
    parser.add_argument('-o', "--output_file", help="Output bam", required=True)

    args = parser.parse_args()

    rewrite_read_groups(args.input_file, args.output_file, args.sample)

If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups.
Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not to be stable across PySAM versions.

I want to assign all read groups to the same sample, where read groups are from different flowcells, but not rewrite all alignments to the same read group, so the Picard AddOrReplaceReadGroups is not the tool I am looking for. I have actually now written something using htsjdk in java to perform the operation I want, but I will leave the question for pysam.

source: Stack Exchange

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

sort: votes


text:
Mini-bulk RNA-seq and scRNA-seq analysis differences?
We have RNA-seq data from libraries prepared using the Smartseq2 single-cell protocol on 500 cells (mini-bulk) / library. The complexity is much better than with single-cells (~14k genes for 1.5M reads).
6 cell types / biological replicates were collected with flow cytometry, and there are 4 control and 5 diseased biological replicates (patients).
We are mostly interested in the differences between control and disease either overall or within each cell type.
Can single-cell tools (e.g. Seurat) be used on this dataset? It is neither bulk nor single-cell data. How would you approach this?

I would process it both ways and check if they are comparable. If not then I would need look into the methods more in detail (To see if they are applicable in this data)

Single-cell tools were developed to address some specific factors unique to scRNA-seq, such as:

a lot of dropouts, so hundreds or thousands of detected genes, not tens of thousands (for mammalian genomes)
hundreds or thousands of replicates
not necessarily well-defined groups

Thus, I would say your experiment is a lot closer to bulk RNA-seq than scRNA-seq.

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


text:
Retrieve detailed gene descriptions
Given a list of gene IDs, how do you retrieve the gene description, summary and other detailed information in R?

MyGene.info
Using the R package for the MyGene.info database:
source("https://bioconductor.org/biocLite.R")
biocLite("mygene")

library(mygene)
gene <- getGene("1017", fields="all")
gene$summary

Multiple genes:
gene <- getGenes(c(1347:1350), fields="summary")

Entrez database
library(devtools)
install_github("ropensci/rentrez")

library(rentrez)
genes <- entrez_summary(db="gene", id=c(718,720))
genes$`720`$summary

Nice question (and answer, I didn't know about MyGene), if you know what you are looking for maybe it is worth to mention biomaRt

Description of how to use biomaRt here. Filter by your list of gene IDs. Get the descriptions etc as attributes. Eg:
> library(biomaRt)
> ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")
> IDs <- c("BRCA2","BRAF")
genedesc <- getBM(attributes=c('external_gene_name','description'), filters = 'external_gene_name', values = IDs, mart =ensembl)
> genedesc
  external_gene_name                                                                      description
1               BRAF B-Raf proto-oncogene, serine/threonine kinase [Source:HGNC Symbol;Acc:HGNC:1097]
2              BRCA2                  BRCA2, DNA repair associated [Source:HGNC Symbol;Acc:HGNC:1101]

Thank you for the answer. I did not find a way to return detailed information  using biomaRt. The description fields are just extended gene names.

The easiest and fastest method (which also doesn’t require internet access) is to use annotation packages, e.g. org.Mm.eg.db:
library(AnnotationDbi)
library(org.Mm.eg.db)

# 1:1 mapping symbol→entrez
entrez_ids <- mapIds(org.Mm.eg.db, gene_symbols, 'ENTREZ', 'SYMBOL')

# 1:n mapping symbol→…
print(columns(org.Mm.eg.db))  # check what you can query
df <- select(org.Mm.eg.db, gene_symbols, c(...), 'SYMBOL')

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

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


text:
Which data format for stranded coverage data
I need tools I am writing to output data in terms of score for each base in the genome. The format needs to be indexed for fast random access.
The obvious choices are bigWig and tabix-index bedgraph (indexed bedGraph seems quicker for computation, bigWig better for visualization), but the problem is that I need the data to be stranded - i.e. I need to record which strand the data belongs to (and perhaps even have a score on the positive and negative strands).
So my choices are:

Invent a new format
Bodge one of the existing formats (this is almost as bad as inventing a new format, if not worse). For example: bedGraph, except with two score columns - one for + and one for -, or have a bedGraph with with a strand column. This doesn't seem ideal because other downstream tools won't be able to process it.
Have seperate files for + and - strand scores. This is unweidly and not great if you sometimes don't care about strand, or want to visualise on a single track.
Use a magic format I've not heard of but that someone is going to point out in the answers...

None of these are ideal. As they either make tool output files incompatible with other tools, or don't encode all the information I want to. Is there a an established format I can use for this?

Since you lose many benefits of creating a new format: I found that having three files — forward, reverse, unstranded/merged — works best.

So far I've been working with three files: a full bed6 with a line per non-zero coverage base. This I index with tabix and use internally. Two bigWigs - forward and reverse, which I use for external programs. The difference in performance is quite impressive. On one benchmark what makes about 18000 queries, querying BAM and processing to stranded depth takes 45 seconds, querying bigWig takes 90 seconds and querying the indexed bed takes 25 seconds.


text:
Count genomic ranges
I have a set of genomic ranges that are potentially overlapping. I want to count the amount of ranges at certain positions using R.
I'm Pretty sure there are good solutions, but I seem to be unable to find them.
Solutions like cut or findIntervals don't achieve what I want as they only count on one vector or accumulate by all values <= break.
Also countMatches {GenomicRanges} doesn't seem to cover it.
Probably one could use Bedtools, but I don't want to leave R.
I could only come up with a hilariously slow solution
# generate test data
testdata <- data.frame(chrom = rep(seq(1,10),10),
                       starts = abs(rnorm(100, mean = 1, sd = 1)) * 1000,
                       ends = abs(rnorm(100, mean = 2, sd = 1)) * 2000)

# make sure that all end coordinates are bigger than start
# this is a requirement of the original data
testdata <- testdata[testdata$ends - testdata$starts > 0,]

# count overlapping ranges on certain positions
count.data <- lapply(unique(testdata$chrom), function(chromosome){
    tmp.inner <- lapply(seq(1,10000, by = 120), function(i){
        sum(testdata$chrom == chromosome & testdata$starts <= i & testdata$ends >= i)
    })
    return(unlist(tmp.inner))
})

# generate a data.frame containing all data
df.count.data <- ldply(count.data, rbind)

# ideally the chromosome will be columns and not rows
t(df.count.data)

So you only want overlaps where the start position is not contained within your interval range? What's the actual biological context to this?

You could start by optimizing the code by using data.table (and probably using vapply instead of lapply). Also if one nucleotide overlaps between ranges then you count that too?

@DevonRyan no what i want is that the respective counting position is with in the range, which is a genomic deletion (hence the biological context)

@Llopis yes i want that counted as well. The list contains overlapping areas and i want to know how many of them are overlapping a certain position. Which is then used later for plotting.

Hi sargas, thanks for your question and welcome to Bioinformatics Stack Exchange. If you have any additional story/context associated with your question (such as a desire to count genomic deletions), it can be helpful to include that in your question. This makes it easier to understand the question, and helps answerers to solve the problem you have rather than the question you're asking (see the XY problem for more information).

GenomicRanges::countOverlaps seems to be what you’re after:
position_range = GRanges(position$chrom, IRanges(position, position, width = 1))
ranges_at_position = countOverlaps(position_ranges, granges)


text:
GOseq analysis with evidence code filter
When using GOseq analysis on RNA-seq data, I often find many 'false positives'. What I mean with that is that some genes in certain categories are not really involved in the process, but are only involved due to 'Electronic Annotation' (IEA). Most of the genes are IEA. A list of all evidence codes is here.
Is it possible to (easily) do GOseq analysis without this evidence group?
I mean, I know how to do it the hard way (make a custom GO-genes annotation in e.g., biomaRt removing all IEA evidence code genes, and perform GOseq with custom genome). But is there also a easy way in GOseq available?

For what it is worth, there have been papers that say that the quality of the IEA is not so bad (as bad as other categories), also the GO says that the evidence codes are not designed to be used as quality filters.

Thanks @Llopis, I mostly agree with their statement about IEA quality. However, I was thinking of ways to prevent getting enriched GO terms such as 'spermatogenesis' in female mice experiments... Another motivation is to get 'real' target genes for further investigation (often you get false positives when further investigating).

Then I recommend to give a look at the topGO package, it should filter most of them out, via taking into account the DAG structure of GO.

@Llopis, do you mean with elim function? That reduces the number of significant GO terms, but won't eliminate false positive genes (targets) for further research.

Yes, I mean the elim and the weight01 methods. I agree that they won't eliminate false positive genes but it will let you focus on genes  on categories that are really relevant. Anyway, let's see if someone came up with a solution, even if it is the long one.

I think the answer to the quesiton is that there is not easy way to do it. But as a comment i'd say that I would never trust a GO analysis to provide "answers", only suggestions for further analysis, and as such any genes in significant categories are always going to have to be examined by a domain expert.

@Llopis, can you direct me to the papers about the quality of IEA and about not using evidence codes for filtering. I am interested in their motivation, thanks.

@b.nota sorry I though I had the references but I don't. I am looking for those articles, if I find them I'll post in a new question and let you know

@b.nota I found the reference! I found that Gene Ontology semantic similarity tools: survey on features and challenges for biological knowledge discovery recommends (under the section "Influence of GO evidence codes") using all the data and evidence codes, because the tools (that make these associations IEA) are getting better to assign a term (and it provides references to semantic similarities with IEA and without IEA). Edit: Now perhaps I should write a question on itself...

@Llopis, thanks for getting back with the reference.

A bit weird to answer my own question, but for who's interested I can show the 'hard' way using biomaRt and then use goseq with gene2cat defined.
So first you need biomaRt to get annotation of the genes. I have ensembl genes from mm10.
library(biomaRt)

table <- read.table(file, header = T, stringsAsFactors = F)

ensembl <- useMart("ENSEMBL_MART_ENSEMBL")

ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)

allGenes <- table$Geneid

GOs <- getBM(attributes=c('ensembl_gene_id', 'go_id', 'name_1006', 'go_linkage_type', 'namespace_1003'),
      filters = 'ensembl_gene_id',
      values = allGenes,
      mart = ensembl)

dim(GOs)
[1] 477082      5

table(GOs$go_linkage_type)
          EXP    IBA     IC    IDA    IEA    IEP    IGI    IMP    IPI    ISA    ISM    ISO    ISS    NAS     ND    RCA    TAS
 36931     36  27918    525  34038 210114    630   6891  30219   6205   4234     20  77365  24673    609  12821     99   3754

# I only want BP and no IEA evidence code
GOs_no_IEA <- GOs[GOs$go_linkage_type != "IEA" & GOs$namespace_1003 == "biological_process",]

dim(GOs_no_IEA)
[1] 111944      5

Then you can continue with goseq.
library(goseq)

All <- read.table(input_All, sep = "\t", header = T)

All_genes <- as.character(All$Geneid)

TopTable <- read.table(input, sep = "\t", header = T)

sig_genes <- as.character(TopTable$Ensembl)

# Make gene vector

gene.vector <- as.integer(All_genes %in% sig_genes)

names(gene.vector) <- All_genes

gene.lengths <- All$Length

pwf <- nullp(gene.vector, "mm10", "ensGene", bias.data = gene.lengths)
head(pwf)

GO.wall <- goseq(pwf, "mm10", "ensGene", gene2cat = GOs_no_IEA[, c(1:2)], use_genes_without_cat = F, test.cats = "GO:BP")


text:
how to merge two fasta based on full fasta header and remove the duplicate
I am wondering if anyone can help to merge two fasta file based on the full header as initially gene ID may be similar in some of the sequences and to remove the duplicate entry based on full fasta header
I do have two multifasta file:
file1:
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

file2:
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

Expected output:
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC

>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

Please [edit] your question and show us i) an example of your input files including "matching" and "non-matching" sequences and ii) the output you would expect from those input files.

Thanks for the edit! Will sequences with the same header be identical? If so, why not extract all unique sequences directly instead of trying to find identical headers?

yes, same header have identical sequences, it will be great to extract unique sequences.

You can use the FastaToTbl and TblToFasta scripts I have posted here for this:
FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa

If you run this on your example, you will get:
$ FastaToTbl file1.fa file2.fa | sort | uniq | TblToFasta > merged.fa
$ cat merged.fa
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATG
GTGCGGGAACGCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGT
TTGGATGTCTCTGGCGACGG
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGC
GAACGCTACGGGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGAT
GTCACCGACGACGGTGGGTG
>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAG
GTCGCGGTCGGCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTC
GACCGCCTCGAAGCCCGCCT
>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAG
GTCGCGGTCGGCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTC
GACCGCCTCGAAGCCAGCGC
>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGC
CAGTTGGTACGGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGAC
GTCGGCCTCGATGTCACGGG

This will remove any cases where both the header and the sequence are identical. If you instead want to do this based on the header only, not taking the sequences into account, you can just process both files and only print sequences from the 2nd if their ID line hasn't been seen before. For example, with awk:
$ awk '{
        if(NR==FNR){
            if(/^>/){seen[$0]++;}
            print
        }
        else{
            if(/^>/){ skip=seen[$0] }
            if(!skip){ print }
        }
       }' file1.fa file2.fa
>AOIT01000069.1:1403-2242 Haloterrigena limicola JCM 13563 contig_69, whole genome shotgun sequence
ATGAGTAACGATACGGTAGCTGGTGACCGCGACCCCGAGGAGGCCCGTGAGATGGTGCGCGAACGCTACG
GGACGATCGCTTCGGACGGTCAGGACTGCTGTGGCGATGTCGGCATCGATGTCACCGACGACGGTGGGTG

>LOPV01000533.1:1144-1992 Haloferax sp. SB29 contig_613, whole genome shotgun sequence
ATGAGTGACGACACCCAGACTGGCACGACAGCCGACGACCGCGATCCAGAGGCGACTCGCCAGTTGGTAC
GGGAACGCTACGCCGGAATCGCATCGAGCGACCAGGGCTGTTGTAGCGACGTCGGCCTCGATGTCACGGG

>AOLJ01000027.1:50569-51417 Haloferax gibbonsii ATCC 33959 contig_27, whole genome shotgun sequence
TCAGTCGTCGAAGGTCGGTTTCCGGGCTTCGATGGTGGCCGACACGAGGTACTCGCCGAGGTCGCGGTCG
GCCTCCCAGTCGCTGATGAACTCAGTACTCTCGTCCTTGGGTGCAATCTCGACCGCCTCGAAGCCAGCGC
>AM774418.1:280304-281149 Halobacterium salinarum R1 plasmid PHS3 complete genome
ATGAGTAATGACAACGAGACGATGGTCGCCGATCGCGATCCCGAGGAGACTCGCGAGATGGTGCGGGAAC
GCTACGCGGGAATCGCGACGAGCGGCCAGGACTGCTGTGGTGACGTCGGTTTGGATGTCTCTGGCGACGG

>AOJK01000067.1:53467-54312 Halorubrum californiensis DSM 19288 contig_67, whole genome shotgun sequence
TCAGTCGTCGCGGGCTGGTTTTCGGGCTTCGATGGTAGCGGAAACGAGGTACTCGCCGAGGTCGCGGTCG
GCGTCCCAGTCGCTGATGAACTCGGTGCTCTCGTCCTTCGGCGCGATCTCGACCGCCTCGAAGCCCGCCT

it will be great to remove identical sequences as well based nucleotide sequence. but it is showing error: $ chmod +x FastaToTbl
$ chmod +x TblToFasta
$FastaToTbl file1.fasta file2.fasta | sort | uniq | TblToFasta > merged.fa
bash: FastaToTbl: command not found...
bash: TblToFasta: command not found...

@Dr.Yogesh yes, that's because you didn't put the scripts in a directory that's in your PATH. If you have them in the current directory, use ./FastaToTbl file1.fa file2.fa | sort | uniq | ./TblToFasta > merged.fa

Thanks, I am trying to copy the script again!

Thanks, it is working perfectly now, it was sensitive to extra space!


text:
Reordering scaffolds according to a reference without a genetic map
I am trying to reorder scaffolds of a rice species, but no genetic map is available right now. Oryza sativa Japonica is a close relative of this rice species. Mummer was used to do a whole genome alignment, and I am trying to reorder scaffolds according the mummer result. But I found it's too complicated to do this by myself, many scaffolds may align to multiple Oryza sativa Japonica chromosomes or multiple loci in the same chromosome.
What software can I use to handle this task?

I've been using ATAC to successfully align different versions of the rice genome assembly. There is a chance it will work for you between two closely related species (potato and tomato work well for example).
It's not massively user friendly, but it's very fast, and it produces long, 'clean' alignments (i.e. almost all 1 to 1 alignments). The output format is a bit cryptic (zero-based start position and length instead of the more usual 1-based start and end positions), but as I said, it should basically output the order for you in about 30 minutes.

That looks like a good solutions. I hopt this will help me to align scafolds

Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.
Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a reference and/or b) rearrange contigs to improve assembly quality.
I use the tutorial as above which uses the GUI for Mauve. I

Could you please expand this and give an example of how Mauve can be used for the OP's purposes?


text:
Does StringTie use exon length for gene abundance or complete gene length?
I'm currently using StringTie with option -A for calculating gene abundance (I'm interested in TPM). While checking the manual, I couldn't find any reference to the record type StringTie uses to calculate this value. Does anyone know if it is total gene length or total effective (i.e. exon) gene length?


text:
No variant found using GATK 4.0 HaplotypeCaller
I am doing variant calling on RNA-seq datasets from wheat which is hexaploid,the binary alignment (BAM) files were created using STAR version 2.6.0c and variant calling was done using GATK 4.0 HaplotypeCaller.The whole pipeline is as follows:
Alignment through STAR 2.6.0c :
STAR --runMode alignReads --runThreadN 20  --outSAMtype BAM  SortedByCoordinate --sjdbOverhang 100  --genomeDir wheat --readFilesIn G1_cleaned_R1.fastq G1_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3

Then AddOrReplaceReadGroups from picard was used with default settings:
 AddOrReplaceReadGroups \
   I=input.bam \
   O=Read_groups_added.bam \
   RGID=4 \
   RGLB=lib1 \
   RGPL=illumina \
   RGPU=unit1 \
   RGSM=20

Finally GATK 4.0 haplotypeCaller was applied:
gatk HaplotypeCaller --native-pair-hmm-threads 28  \
    -I Read_groups_added.bam -R wheat_ref.fa -O final.vcf

The command ran successfully as shown below:
20:33:01.018 INFO  ProgressMeter - TGACv1_scaffold_470693_5DS:301            102.2              45108670         441320.0
20:33:11.039 INFO  ProgressMeter - TGACv1_scaffold_724195_U:301            102.4              45111270         440625.5
20:33:21.072 INFO  ProgressMeter - TGACv1_scaffold_724579_U:301            102.5              45113880         439932.4
20:33:31.108 INFO  ProgressMeter - TGACv1_scaffold_248869_3B:301            102.7              45116440         439240.9
20:33:41.143 INFO  ProgressMeter - TGACv1_scaffold_248918_3B:301            102.9              45118990         438551.6
20:33:51.150 INFO  ProgressMeter - TGACv1_scaffold_306137_4AL:301            103.0              45121550         437866.7
20:34:01.176 INFO  ProgressMeter - TGACv1_scaffold_193559_2DS:301            103.2              45124100         437182.5
20:34:03.134 INFO  HaplotypeCaller - 100112154 read(s) filtered by: ((((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter) AND WellformedReadFilter)
100112154 read(s) filtered by: (((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter) AND GoodCigarReadFilter)
100112154 read(s) filtered by: ((((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter) AND NonZeroReferenceLengthAlignmentReadFilter)
100112154 read(s) filtered by: (((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter) AND PassesVendorQualityCheckReadFilter)
100112154 read(s) filtered by: ((((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter) AND NotDuplicateReadFilter)
100112154 read(s) filtered by: (((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter) AND NotSecondaryAlignmentReadFilter)
100112154 read(s) filtered by: ((MappingQualityReadFilter AND MappingQualityAvailableReadFilter) AND MappedReadFilter)
100112154 read(s) filtered by: (MappingQualityReadFilter AND MappingQualityAvailableReadFilter)
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

20:34:03.134 INFO  ProgressMeter - TGACv1_scaffold_725922_U:301            103.2              45124592         437049.1
20:34:03.134 INFO  ProgressMeter - Traversal complete. Processed 45124592 total regions in 103.2 minutes.
20:34:04.749 INFO  VectorLoglessPairHMM - Time spent in setup for JNI call : 0.0
20:34:04.749 INFO  PairHMM - Total compute time in PairHMM computeLogLikelihoods() : 0.0
20:34:04.749 INFO  SmithWatermanAligner - Total compute time in java Smith-Waterman : 0.00 sec
20:34:04.749 INFO  HaplotypeCaller - Shutting down engine
[July 3, 2018 8:34:04 PM PKT]   org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller done. Elapsed time: 103.99 minutes.
Runtime.totalMemory()=4789370880

but still the generated vcf file i.e final.vcf does not contain any variants as shown below:
[dellemc@localhost bilal]$ tail final.vcf
##contig=<ID=TGACv1_scaffold_725916_U,length=500>
##contig=<ID=TGACv1_scaffold_725917_U,length=500>
##contig=<ID=TGACv1_scaffold_725918_U,length=500>
##contig=<ID=TGACv1_scaffold_725919_U,length=500>
##contig=<ID=TGACv1_scaffold_725920_U,length=500>
##contig=<ID=TGACv1_scaffold_725921_U,length=500>
##contig=<ID=TGACv1_scaffold_725922_U,length=500>
##contig=<ID=TGACv1_scaffold_731344_7BL,length=500>
##source=HaplotypeCaller
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  20

Is there any mistake in the haplotypeCaller command or any argument is missing ?

All those compute times of 0.0 seem a little suspicious.

I have edited the question to add more details.

@swbarnes could you tell about the type of suspicions?

There are steps of the 'best practices' workflow missing, which I think is leading to the majority of reads being filtered by the haplotype caller, notice these lines in the output:
21066679 read(s) filtered by: MappingQualityReadFilter
79045475 read(s) filtered by: MappingQualityAvailableReadFilter

That's >100m reads being filtered by the HC.
Specifically, you're missing SplitNCigarReads, which separates reads that bridge a splice junction (padded by Ns in the STAR alignment) and reassigns the mapping quality. At the moment this element of the workflow is not validated or available in GATK4 (RNA-Seq variant calling is largely abandonware in newer versions of GATK, AFAIK), but in the best practices WDL utilises GATK3.5 - see https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/master/rna-germline-variant-calling.wdl, line 440-477.
There are other elements of the recommended workflow missing too, but I think this is the critical bit. I won't harp on here about the unsuitability of RNA-Seq data for variant calling, as I'm all too aware there often isn't much option.

What is meant by the option --limitOutSJcollapsed in the STAR command given in the link you suggested https://github.com/gatk-workflows/gatk3-4-rnaseq-germline-snps-indels/blob/a7a890b38af1309edc040ac89814df37a57f77b7/rna-germline-variant-calling.wdl#L345

A 'collapsed' junction is what gets written to the SJ.out.tab file during a run of STAR (the collapsing here refers to the removal of duplicate junctions). This option merely limits the maximum number of these junctions that are written. It defaults to 1,000,000.

As well as the issues @sjcockell mentions, you have another potential problem you've got a hexaploid organism, by default the HaplotypeCaller (HC) is hardwired to a ploidy of 2, you need to set this to 6 to match your organisms ploidy using -ploidy 6 as an argument.  Since the ploidy is a key term in the underlying Bayesian stats the HC employs (see this link for an explanation) it's likely that some of the lower allelic frequencies that hexaploidy would generate will mean quite a few variants won't be called as they'd be statistically unlikely for a diploid organism.

You might also want to consider the GATK support forum too, as they are very helpful on GATK matters: https://gatkforums.broadinstitute.org/gatk/categories/gatk-support-forum


text:
Cannot install chromosomer
I am trying to install chromosomer but I fail. Can anybody help me, please?
$ pip install chromosomer

Collecting chromosomer
  Could not find a version that satisfies the requirement chromosomer (from versions: )
No matching distribution found for chromosomer

The output python:
$ python

Python 3.6.0 |Anaconda custom (x86_64)| (default, Dec 23 2016, 13:19:00)
[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] on darwin
Type "help", "copyright", "credits" or "license" for more information.

The command line output from pip isn’t relevant for the question, I’ve thus removed. pip --version would have been more relevant. Also, be sure to read the formatting help for questions, especially the part about preformatted text (which applies to console output).

Chromosomer only exists for Python 2. You should thus be able to install it via
pip2 install chromosomer

But:
Your Python installation is a bit screwed up: there are Python 2 and Python 3, which are unfortunately incompatible.
On your system python and pip seem to be aliases for Python 3, which I’d strongly recommend against (lots of tools will break). Instead, they should alias Python 2. Python 3 should be invoked via python3/pip3.

That's a pretty sweeping statement. The command python defaults to python3 on my two Arch systems and various Ubuntu boxes we use at work and I've never had issues with it. What makes you say that using python2 as the default is better? And on what systems?

@terdon It's from a PEP recommendation. And failing to do so, while not catastrophic, definitely breaks many scripts that assume python = Python 2. Not a fan of Arch's/Ubuntu's decision to ignore that unnecessarily.

Well, I'd argue that using a shebang line calling python instead of python[23] is a very bad design decision since the python devs, in their infinite wisdom, made the two versions incompatible. Either way, your advice is dangerous since, presumably, the OP has a working system with python3 as the default and changing that will break lots of tools while leaving it as is will not.

@terdon You realise that there are Python scripts that predate the existence of Python 3, right? The all will, by necessity, have python in their shebang line. And fixing the system shouldn't break anything but, yes, it will require reinstalling all modules.


text:
Human body sites knowledge base
I'm looking for a lightweight knowledgebase describing the human body to annotate disease sites.
I do not need a great level of detail, I just need a kind of basic organs/sub-organs taxonomy. I checked out the Foundational Model of Anatomy ontology, but it's over-kill for my purposes and really cumbersome (Protege fails when it tries to load FMA.owl).
Is there a simple online knowlegebase available for this?

I did manage to open FMA.owl in Protégé, by the way. It took a while and a considerable amount of CPU and RAM but it did work eventually.

Ok! Probably I have no enough HW resources to load FMA and then extracting a small subset from it. That's the reason why I was looking for something smaller.

I suggest you post a new question asking how to extract a subset of an existing ontology. Give us a link to FMA, explain what you need (for instance, how deep you want to go in the ontology's branches) and ask for help making a "slim" ontology from it.

I am the developer of Uberon and I would be happy to help you with what you need, either from Uberon, or from the FMA.
You mentioned you need something simpler than FMA. There are a variety of tools for creating custom subsets. Additionally, some ontologies provide ready-made subsets for particular purposes.
For example, we make a subset called 'basic.obo' from Uberon:
http://uberon.github.io/downloads.html#subsets
You can examine this directly in either oboedit or Protege5 (note Protege5 parses obo so you don't need to convert).
From your question, you didn't indicate you needed programmatic access, but if you do there are a variety of options, as well as more programmer-friendly JSON exports.
Most ontologies have a number of different options for web-browsing, you can see some of the options for Uberon here:
http://obofoundry.org/ontology/uberon.html
You mentioned annotating disease sites. There may be a tools (either desktop or web-based) for helping with this - e.g. providing autocomplete over a desired subset of an ontology.
You may also want to ask your question on a mailing list like obo-discuss. Unfortunately many ontology developers don't read stackexchange, but most are always keen to hear of community requirements, to help them provide better and more useful products.

Search in bioportal: For example:  SNOMED, UBERON.
Please, note that FMA is also available as RDF, csv , etc...
You could also try to convert FMA to obo and try to open it with OBO-edit

SNOMEDCT seems to be a full ontology for any type of clinical data, not only anatomical (e.g. it includes terms such as "smoker" or "observation") which would make it far more complex than the FMA the OP tried. Likewise, Oberon seems to be a cross-species ontology which, again would be more, not less complex than FMA. If I understand correctly, the OP wants something less complete and less complex than FMA.

Yes @terdon, you are right. I'm looking for something simpler and smaller than FMA.

@Pierre, I'm afraid that Oboedit performs even worse than Protege


text:
Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header?
Is it possible to find out the PacBio chemistry from the PacBio Sequel BAM's header? E.g. whether it was produced with P5-C3 or P6-C4? I checked mine but I haven't found it but maybe it's hidden in one of the tags whose meaning I don't know yet, like SEQUENCINGKIT=100-861-800?

According to PacBio Sequel BAM format manual, you can find information on the chemistry version in the header entries @RG (read group), tag DS (description), entry BASECALLERVERSION .
In case of future very PacBio-specific questions, you can contact Liz Tseng on Twitter (@magdoll), one of the most responsive PacBio developers. They also have a dedicated SeqAnswers forum, which is directly monitored by them and you are very likely to get a fast and accurate answer.

I have 5.0.0.6236 for BASECALLERVERSION. Does it mean chemistry P5-C3 was used?

This surely means chemistry is higher than P4-C2 / 2.0, which is not supported in pbbam. The logical answer would be then yes, P5-C3 chemistry (or higher) was used.


text:
Pooling data in metagenome assembly
I have 12 human gut microbiome WGS Nextseq reads (151 bp paired end). What will be an effective strategy to assemble a metagenome?
Let us say I have already filtered the fastq for quality, adapter sequence and host contamination (human, in this case).
1) Should I concatenate all the R1 reads as one single R1 read and one Single R2 read?
cat Sample[1..12].R1 > Single_R1.fastq
cat Sample[1..12].R2 > Single_R2.fastq

and then use Diginorm to normalize the Single_R1.fastq and Single_R2.fastq. Subsequently, feed these fastq files into any metagenome assembler such as Megahit, MetaSPAdes?
Normalize the output by using CD-HIT or similar tool to remove duplicates and filter by contig length.
OR
2) Perform metagenome assembly for each of the samples individually after applying filtering, removing adapters and host contamination.
R1=(*_R1_001.filtered.fastq)
R2=(*_R2_001.filtered.fastq)
for ((i=0;i<=${#R1[@]};i++)); do
  /bin/metagenome-assembler -1 "${R1[i]}" -2 "${R2[i]}" -o ${R1[i]%.*}.contigs.fa;
done

Followed by combining all the contigs.fa into one mega_contigs.fa
cat *.contigs.fa > Mega_contigs.fa

and use CD-HIT or similar tool to remove duplicates.

I think assembling data separately is better idea, since there will more noise in combined signals of different individuals and also resource wise. However, I guess you should not just cat assemblies together, you want to avoid redundancy (discard what is shared between individuals), not entirely sure how to achieve that though.

With metagenomes, there's several strategies and no one-size-fits-all.
 You can spend a lot of time playing with the data to try to get the "best assembly."  But it depends on what you're trying to do, e.g. whether or not you are mostly interested in the highest abundant species or are interested in the less abundant ones too.
But to propose a strategy that seems reasonable and hasn't been mentioned yet, my initial thoughts are that since the gut microbiome isn't complex, you can likely get decent assemblies from each run separately.  Then you can run a contig clustering algorithm (e.g. MetaBAT) to bin the contigs into genomes.  Next compare the bins from the various assemblies (i.e. cluster the clusters).
It may or may not be worth the trouble to try to combine the related bins' contigs into a consensus (pan)genome sequence.  If you want to, methods that have been used successfully include shredding the contigs into overlapping chunks and assembling with an OLC assembler.  But if you're interested in KO:read-counts (for example), a consensus genome sequence isn't really required.  Cross-mapping the unassembled reads can be useful too (i.e. a sample may have a low-abundance OTU that doesn't produce good contigs, but it can be associated with an OTU that is better assembled in another sample).
If you provide more detail on your goal, I can comment further.

If your goal is to bin the resulting contigs into genomes, then you should do option #1, pooling the reads and assembling into one set of contigs.

I also think pooling is the better option, followed by partitioning by coverage/taxonomy of contigs.
Maybe check out BlobTools, which helps with filtering read-pairs by taxonomy of contigs they contribute to and also does nice visualisations of assemblies.
Workflow B seems to be what you want.
Disclaimer: I am the developer of this tool.


text:
5'UTR and 3'UTR annotation in yeast
I am working on a project in which I need to compute several parameters (such GC content and length) of 5'UTR and 3'UTR sequences of Saccharomyces cerevisiae yeast genes.
The problem is finding a proper annotation for these regions in yeast. I have tried with BiomaRt and these sequences are not available for yeast. I have also tried with UCSC table browser to obtain a BED file with 5'UTR or 3'UTR sequences but it is returning only coordinates of tRNAs, pseudogenes, ncRNAs.
In summary, is there a straight-forward way to obtain an annotation for 5'UTR and 3'UTR regions in Saccharomyces cerevisiae yeast?

Sorry, is for Saccharomyces cerevisiae. I am going to edit the question just in case.

isn't there an equivalent of AceView for Yeast ? https://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/index.html

I am unaware of any "official" or gold-standard UTR annotations in S. cerevisiae.
One option is to use the annotations from the TIF-Seq publication (Pelechano et al. 2013).
The GSE39128_tsedall.txt.gz file contains the major isoforms identified. It would be up to you to computational associate each transcript with a given gene. It is also up to you to determine which major isoforms are present in your dataset. Furthermore based on my own experience there is often a significant number of discrepancies between what I see in an NGS dataset (e.g. RNA-Seq) and what is called as a possible isoform in these annotations. However, a lot of this could be related to the inherent technological limitations of the method.
If you look at the SGD genome browser they provide some additional UTR annotations. Derived from Nagalakshmi et al. 2008.
You can find the relevant annotations for those UTR's in their supplemental data section. The file is 1158441_tables_s2_to_s6.zip and inside you will find TableS4 which seems to contain what you seek.
There might be more resources but these are the two main ones I'm aware of.

I found the following two files in https://downloads.yeastgenome.org/sequence/S288C_reference/:

SGD_all_ORFs_3prime_UTRs.fsa
SGD_all_ORFs_5prime_UTRs.fsa

According to the README files in the same directory, these are (the README for the 5' file is equivalent):

SGD_all_ORFs_3prime_UTRs.README
Information about the SGD_all_ORFs_3prime_UTRs.fsa file.
This file contains sequences for predicted 3' UTR sequences for all
  characterized, uncharacterized and dubious ORFs in S288c.
3' UTR coordinates were retrieved from YeastMine
  (yeastmine.yeastgenome.org) on July 3, 2017 using the stop coordinates
  of both the predicted transcript and gene sequences. A BED file with
  these coordinates was made. BED files were fixed on Sept 22, 2017 to
  account for strandedness. Sequences were batch downloaded again in
  FASTA format from USCS Genome Table Browser on Sept 22, 2017
  (https://genome.ucsc.edu/cgi-bin/hgTables?hgsid=604915183_aUc9FZ1vwyAccGti1uI0wteYEM7x)
File created: 2017-08-25

Those therefore seem to be the 3' and 5' UTR sequences of all yeast genes in FASTA format and should be exactly what you need.

The link is broken

@Pengyao yes, apparently yeastdownloads.org are using a bad SSL certificate. I guess they'll fix it eventually.


text:
Calculating the charge of a peptide computationally
I was wondering how I can calculate the charge of a protein peptide (e.g. "RKTTLVPNTQTASPR") computationally in R or another tool.

Have you done any search already? What have you found?

A quick google search turns up protcalc, which is able to give a nice pH-dependent table of peptide charges (yours ranges from 3.1 at pH 4 to 1.5 at pH 10). It's on sourceforge, so hopefully the source code (and maybe a publication) is somewhere in there. Granted, this isn't R, but it's not clear from your post how necessary that really is.
BTW, in general "best and most accurate" turns out to be highly subjective (at least the "best" part of that).

Actually somewhere I read that a method to calculate the charge was taking in mind the pkAs of the free amino acids and not the residues on the chain. That's why I wrote "best". I think that the pkA change when the a.a is on the chain. Maybe I should have mentioned it. Also R environment will be preferred.

I would be surprised if pH still didn't affect things, though I agree that using a calculation based on free amino acid charges is going to yield wrong values.

Taking in mind the pH is standard but pH value is known and can be found. The thing is if there is an algorithm to calculate the pkA of the non free a.a

@J.Doe asking for "the best" is not a good question for this site since "the best" is almost always going to be a matter of opinion, and/or depend on the specifics of the sequence you are analyzing. See https://bioinformatics.stackexchange.com/help/dont-ask

Another quick Google search reveals the Peptides R package.
library(Peptides)
s <- "RKTTLVPNTQTASPR"
charge(s)

[1] 2.997683

Optional arguments to charge() are pH (default = 7) and pKscale (choice of nine, default = "Lehninger"). See ?charge for details.
Another (non-R) option is the EMBOSS suite program iep. I highly recommend installing EMBOSS, it provides a suite of useful command-line tools with output that can be easily read and processed by other sofware, including R.

I've used PropKa in the past to get isoelectric points. Pretty simple to use:
https://github.com/jensengroup/propka-3.1


text:
Find most similar genes to a set of genes of interest in time series RNA-seq data
I have gene expression data (RNA-seq) for 30 different time points (from 0 to 60 min each 2 min).
I have a set of 8 genes that behave similarly (although not identically) and I want to find the top X genes that behave most similarly to this set of genes in terms of their time-series expression profile.
I have tried so far to do unsupervised clustering (mfuzz) and check for genes that lie in the same cluster as my genes of interest. However, this approach relies on specifying a predefined number of clusters and since these genes are not exactly equal, sometimes they don't fall in the same cluster, depending on initial parameters.
Because of this, I have thought about doing it the other way around and look for genes similar to my set of genes. However, I haven't found any method/resource for doing that so far.
So, is there a way to find genes with a time-series profile similar to a predefined set of genes?

I'm not an expert but hierarchical clustering could be an alternative start

DPGP (McDowell et al, 2018) uses a non-parametric clustering method, so should find clusters without you needing to specify cluster size or number of clusters.  It is more or less explicitly designed for clustering timeseries gene expression data to find clusters of co-expressed genes.
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005896
SplineCluster (Heard et al. 2006) uses a different method, but for a similar purpose and in a similarly non-parametric way, and worked well when we used it a few years ago with microarray data.

I already have a tool for doing the clustering (mFuzz). What I am looking for is to take 8 genes and look for the top X most similar genes to those genes. But I will try this DPGP, seems pretty good for the data I have!

You can use GSVA to summarize the information of your 8 genes and look for the time-serie change using limma or DESeq2.
If you have enough samples for each time point you could also do a weighted gene co-expression network analysis (WGCNA), and see if those clusters in each time point are co-expressed similarly (But it can be tricky and you will need to define a lower bound of the clusters to 8 or so). You could also divide your time series in two point at the beginning and at the end, and see if there are differences or not.
If you want to get a list of genes similar to your set of genes you can use GOSemSim or BioCor (Disclaimer I am a contributor to the first, and author of the second), which can calculate similarities between genes based on the gene ontology and on pathways respectively. For this approach you'll need to calculate the similarity between all the genes and then select those that are closer to your genes of interest.
I'm not aware of any database/website with time-series profiles of genes (I assume that if there would be any would be on humans or some simple model organism).


text:
tblastn not providing full alignment even on exact-matching sequences
I am using tblastn to search a polypeptide seq against transcripts.  One of those transcripts contains the residues for that exact polypeptide once translated, yet BLAST only gets 88% of it rather than it being full length.
Here's the peptide:
>Syntaxin
MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDD
KMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQID
YRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRD
LHDMFMDMAMLVESQGEMIDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGV
IGGTLG

And the transcript.
Here's the first-frame translation of the transcript which contains the exact polypeptide above:
frame +1
> - 669
PATRNITIVSIYDYH*NNDLKCVFTE*MEILNRIWILLNSHSP*AQRLHAI**MTKDRLAALKAAQSDDDD
NDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDDKMKEELEELMSEIKKNAN
KVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGK
TTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEM
IDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG***PSSLPLVYS
LIQLFYAACTTAASNIWTHHNRNSYYYNSNSNSNSNCYSNYFSNYYSNRKHYYYYHYYYYYYYYYYYYYYY
YYYYYYDYYYYYYYYYYYYYYYYYCYYYYYYYYYHFHYYYHHYYYHYYYYHYYYHYHHHHHYHYHYHYHHH
HHYHYHYHYHYHYYYYYYYYYYYYYYYYYYYYYYCYYYYYYTFSLVISLYGSSEM*MRLLVLSLYTFYGLF
YPLQQKGSFFF*SLQQTGTGFYNFFSYCQLHNLSIPSPSYFF*EWVFSFLKKHLIEFFCFLILLKFVLSFS
PIISLIISVCRKKISRIFYTIYLKMS**NHYLLF

Yet when I do a local tblastn (using the modern blast+) or use NCBI's pairwise blast tool it doesn't return a full-length match.
Why does it leave off the last 12% of the sequence match?
Here's the alignment:

That's because by default, blast masks low complexity regions. With the online blast2seqs tool, you can disable this by clicking on the "Algorithm options" button:

And then deselect the "Low complexity regions" filter:

This will let you align the entire sequence as expected:
Alignment statistics for match #1
Score   Expect  Method  Identities  Positives   Gaps
Frame
588 bits(1516)  0.0 Compositional matrix adjust.    290/290(100%)   290/290(100%)   0/290(0%)
+1
Query  1     MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  60
             MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH
Sbjct  154   MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  333

Query  61    SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  120
             SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT
Sbjct  334   SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  513

Query  121   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  180
             QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT
Sbjct  514   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  693

Query  181   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  240
             QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE
Sbjct  694   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  873

Query  241   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  290
             QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG
Sbjct  874   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  1023

This is the part of your transcript sequence that was masked:
acgctaggctgaccctcctccctccctctagtctactccttaatacaattgttctacgccgcgtgcacc

I am guessing it is the repeated ccct nucleotides that resulted in the sequence being masked, but I am not sure.


text:
What database do gene IDs starting with sp_v2 refer to?
I am a student trying to analyze GEO2r datas for one of my courses.
The IDs given in the output are different for different series. I need to convert all of them to a similar format.
In this process I encountered the following type of ID which I don't know its origin:
”ID"    "logFC"
"SP_v2 4634"    "-0.9897758"
"SP_v2 3382"    "-0.8391782"
"SP_v2 4210"    "-1.1693583"
"SP_v2 2117"    "-1.0504727"
"SP_v2 3488"    "-0.9756444"
"SP_v2 1128"    "-0.8289103"
"SP_v2 2735"    "-0.8629999"
...

Each one of the rows represent a single gene.
My question is that what is this ID? The GEO accession is GSE97750.

Looks like some home made accession code, from which GEO set is it?

The GEO accession is GSE97750

Yeah it is a custom array. Annotation can be found in the platform (GPL22166) data table. See answer by sjockell.

This is possibly a wild goose chase, but a lot of searching led me to this sample:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2536852
Which appears to exhibit accessions of the right format. This in turn leads to the platform:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL22166
Which is a custom spotted cDNA array, which the metadata seems to suggest is a human platform.
The full annotation table gives gene symbols and Entrez IDs for most of the probes on the array:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL22166&id=53508&db=GeoDb_blob144
To take the example of the first line of your results:
ID          ORF     Entrez gene SEQUENCE
SP_v2 4634  ETNK1   55500       AAAGCAGCTTCATCTTTCAAAATTGATTTGCTCTGGTTTT

The Entrez Gene record for that ID matches up:
https://www.ncbi.nlm.nih.gov/gene/?term=55500


text:
Does Prokka or gff3 change the sequence / subject ID of fasta files?
I used prokka to create a database for some fasta files but I noticed a strange difference between the prokka and fasta files. Normally the seqid for the gff3 output for prokka has seqid's that match the sequence id's of the fasta files. The files I just produced with prokka have a different seqid than my fasta files. My fasta files look like this
>NODE_108_length_645_cov_0.679537_ID_215

However the prokka output gff3 file for the same fasta file uses a seqid like this
gnl|X|JHIOGLGG_28

What would cause this difference? And do the last digits still correlate to the same seqid's even though the string before it has changed?

Prokka is an annotation tool. It takes your contigs file with fasta headers >NODE_108...and searchs for ORFs and their putative functions. Prokka produces submission ready output, the new fasta headers you see >JHIOGLGG_28...are locus tags and gene description from the annotation and are not related to contig node data. Please see the github information page for more info.


text:
what is the difference between local identity and global identity for homolog finding?
We believe that if after running blast, the global identity between a resulting sequence and our query is at least 30%, we can say that two sequences are homologs. what is the difference between local identity and global identity and how can we calculate them?
the following file is the global alignment result of one of my PSI-BLAST hits against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 264
# Identity:      85/264 (32.2%)
# Similarity:   125/264 (47.3%)
# Gaps:          29/264 (11.0%)
# Score: 376.5
#
#
#=======================================

NP_418577.1        1 -----------MAEMKNLKIE-------VVRYNPEVDTAPHSAFYEVPYD     32
                                ..|..:..|:       :.|::||||..|....|:|...
WP_036312822.      1 MTIVDSGAPADTQEANDSGIQSYLVTFIIRRFDPEVDAEPRWVDYDVEMY     50

NP_418577.1       33 ATTSLLDALGYIKDNLAPDLSYRWSCRMAICGSCGMMVNNVPKLACKTFL     82
                     .|..:||||..||.::...||:|.||...||||..|.:|...:|||||.:
WP_036312822.     51 PTDRVLDALHRIKWDVDGTLSFRRSCAHGICGSDAMRINGRNRLACKTLI    100

NP_418577.1       83 R--DYTDGMKVEALANFPIERDLVVDMTHFIESLEAIKPYIIGNSRTADQ    130
                     :  |.:..:.|||:...|:|:||:|||..|.||...::|::...|.....
WP_036312822.    101 KDLDISKPIYVEAIKGLPLEKDLIVDMDPFFESFRDVQPFLQPKSAPEPG    150

NP_418577.1      131 GTNIQTPAQMAKYHQFSGCINCGLCYAACPQFGLNPEFIGPAAITLAHRY    180
                     ....|:....|.|...:.||.|..|.::||.|..:.::.|||||..|||:
WP_036312822.    151 KERFQSIKDRAVYDDTTKCILCAACTSSCPVFWTDGQYFGPAAIVNAHRF    200

NP_418577.1      181 NEDSRDHGKKERMAQLNSQNGVWSCTFVGYCSEVCPKHVDPAAAIQQGKV    230
                     ..||||.....|:..||.:.|||.|.....|:|.||:.::...||.:.|.
WP_036312822.    201 IFDSRDDAADVRLDILNDKEGVWRCRTTFNCTEACPRGIEITKAIAEVKQ    250

NP_418577.1      231 ESSKDFLIATLKPR    244
                     ...:.
WP_036312822.    251 AVLRG---------    255

#---------------------------------------
#---------------------------------------
the following file is the local alignment of that hit against query:
Aligned_sequences: 2
# 1: NP_418577.1
# 2: WP_036312822.1
# Matrix: EBLOSUM62
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 219
# Identity:      83/219 (37.9%)
# Similarity:   119/219 (54.3%)
# Gaps:           2/219 ( 0.9%)
# Score: 390.5
#
#
#=======================================

NP_418577.1       13 RYNPEVDTAPHSAFYEVPYDATTSLLDALGYIKDNLAPDLSYRWSCRMAI     62
                     |::||||..|....|:|....|..:||||..||.::...||:|.||...|
WP_036312822.     31 RFDPEVDAEPRWVDYDVEMYPTDRVLDALHRIKWDVDGTLSFRRSCAHGI     80

NP_418577.1       63 CGSCGMMVNNVPKLACKTFLR--DYTDGMKVEALANFPIERDLVVDMTHF    110
                     |||..|.:|...:|||||.::  |.:..:.|||:...|:|:||:|||..|
WP_036312822.     81 CGSDAMRINGRNRLACKTLIKDLDISKPIYVEAIKGLPLEKDLIVDMDPF    130

NP_418577.1      111 IESLEAIKPYIIGNSRTADQGTNIQTPAQMAKYHQFSGCINCGLCYAACP    160
                     .||...::|::...|.........|:....|.|...:.||.|..|.::||
WP_036312822.    131 FESFRDVQPFLQPKSAPEPGKERFQSIKDRAVYDDTTKCILCAACTSSCP    180

NP_418577.1      161 QFGLNPEFIGPAAITLAHRYNEDSRDHGKKERMAQLNSQNGVWSCTFVGY    210
                     .|..:.::.|||||..|||:..||||.....|:..||.:.|||.|.....
WP_036312822.    181 VFWTDGQYFGPAAIVNAHRFIFDSRDDAADVRLDILNDKEGVWRCRTTFN    230

NP_418577.1      211 CSEVCPKHVDPAAAIQQGK    229
                     |:|.||:.::...||.:.|
WP_036312822.    231 CTEACPRGIEITKAIAEVK    249

#---------------------------------------
#---------------------------------------

BLAST (Basic Local Alignment Search Tool) doesn't give you global alignments, so you may need to reconsider what you're doing/asking.
Furthermore, "homolog" has a very specific definition. Homologs are genes that are both descended from the same ancestral sequence. It doesn't, as many often mistake, mean 'similarity'. 30% sequence identity, does not imply homology by itself. You need more information.

Take a look here (http://www.ebi.ac.uk/Tools/psa/) if you want to find tools for local and global alignment.

no it's not my problem. Maybe I couldn't describe it well. I have run PSI-BLAST for my sequence and I 've got some sequences that I believe they are homologs of my sequence, but I'm not sure if they are really homologs or not and which one of these sequences are real homologs. For making sure about this homology, I must check local and global identity and I had done it. Maybe a part of a sequence with mine, has a good identity, but I can't say these two sequences are homologs. They must have globally good identity

Perhaps I'm wrong but "They must have globally good identity" still isn't really enough to say that these sequences are homologs. You need ancestry information for the organisms you're interested in.
If you really want to know orthology/homology, I'd suggest using a dedicated orthology finding program like OrthAgogue, OrthoMCL, or better yet ROARY if your sequences are bacterial.

No they don't. It all depends on how closely related the species are. Please [edit] your question and add more information. There is no magix sequence simlarity threshold that defines homology.

how much the species should be close?

There is no magic formula. For closely related species I would expect a sequence identity of far more than 30% for a true homolog but for distant species, I could also expect far lower. As I told you in my answer to your previous question, it is far more important to look at the conservation of specific residues than to look at the overall sequence similarity value.

Please don't post images of text. They just make the question slower to load, they cannot be edited and cannot be searched. Instead, just copy/paste the text directly into your question and use the formatting tools to make it appear like code. Also, the alignment is useless if we don't i) know the species being aligned and ii) know the e-value of the alignment.

If you want a purely alignment based approach to begin to infer homology (note I say infer, as we stated in the comments, this isn't sufficient to say 2 sequences are homologs), I would use a Multiple Sequence Alignment (MSA) instead of a pairwise like you've done so far, and include as many similar proteins as you can from BLAST.
The choice of algorithm for the MSA probably won't matter too much, as you will be looking for conserved regions anyway.
Once you have your MSA, you can easily spot conserved domains. If anything is known about this protein structurally, it's quite trivial to work out whether a conserved region is a functional domain (active site, that kind of thing).
With this information, you would be somewhat justified in calling these proteins structural homologs at least if they share domains, but, to reiterate, without doing phylogenetic analyses, you cannot call them true homologs.

Defining homology isn't as simple as an identity cutoff, however, you're right that homology seems likely given a certain identity. This is where e-value comes into play. See:
https://en.wikipedia.org/wiki/Expected_value
To answer the question, global-global requires that the alignment account for every base in the query and the target sequence. local-global requires a the full length of the query on a part of the target. local-local allows only the highest scoring part of the alignment between the query and target to be reported. i.e. one functional domain against another. these hits are often called HSPs.


text:
Obtaining all protein sequences with a particular domain architecture from Pfam
I want to get the alignment of chain A of 1kf6 (PDB ID) from the pfam database here. This protein chain has two main domains (FAD_binding_2 and    Succ_DH_flav_C). In pfam there is a link to one of these domains and after clicking in one of the mentioned domains in the table below the page, another page comes at the top of which there is a link to architectures.
For example, if we click on "FAD_binding_2" in the table, there are 220 architectures that have this domain and after clicking on that button (at the top of page), a page opens here. In this page 7471 sequences with the following architecture: FAD_binding_2, Succ_DH_flav_C, resembles most to my chain (because it has both FAD_binding_2 and Succ_DH_flav_C domains).
If we want to see all 7471 sequences with that architecture we can click on the show button under that architecture and finally we can see all of that 7471 architecture resembling our chain A sequence. Down the page we can see other architectures that all have FAD_binding_2, but some of them do not have both FAD_binding_2 and Succ_DH_flav_C domains (or some of them have other extra domains that we are not interested in) and so we do not need their sequences and their alignments.
When we want to get the alignment of full sequences at the left side of the page, we can see that we can just get the alignment of all 15696 sequences and it seems that there is no way to get the alignment of our favourite 7471 sequences. I would like to know if there is any way to get the alignment or at least the fasta file of our favourite 7471 sequences (which have both of our favourite domains not just one of them)?

I believe you can do this using SQL files of the database.
First you need to get the integer accession of this architecture. architecture table has the following columns:

auto_architecture (integer id, this is what we need)
architecture (a string describing the architecture)
type_example (example sequence)
no_seqs (number of sequences)
architecture_acc (accessions of the architectures)

You need the following line:
3900719357      FAD_binding_2~Succ_DH_flav_C    Z9JRB3  7471    PF00890 PF02910

Now you need to get the sequences with this architecture id. First, download the file with all the sequences (warning, this one is more than 7Gb). With this one you need column #1 (pfam sequence accesion), #12 (sequence itself) and #16 (architecture integer id).
If you do not want to import files in the database, you can use awk to export all the sequences to fasta format.
zcat pfamseq.txt.gz | awk -F\\t '{ if ($16=="3900719357") print ">" $1 "\n" $12}' > sequences.fasta

Of couse, if you need to perform multiple queries like this, it is much easier to import tables into the database. Then you can perform both queries in a single SQL operation (assuming you know the architecture string):
SELECT pfamseq_acc, sequence
    FROM architecture
    JOIN pfamseq ON pfamseq.auto_architecture=architecture.auto_architecture
    WHERE architecture="FAD_binding_2~Succ_DH_flav_C";

(I haven't tried this SQL query, so maybe you need to correct it a bit)

I think the idea of not importing 7Gb database is better because I have a lot of protein sequences that I want to do similar task on. I prefer trying the second idea, but I don't know why I couldn't get the integer accession of this architecture from above "architecture table" link

Not sure I understand your comment well. In my example I assume you only have the architecture string. If you already have your auto_architecture accession, you can simply write SELECT pfamseq_acc, sequence FROM pfamseq WHERE auto_architecture=3900719357.

An alternative way is using the Hmmer website, and do a hmmsearch with the Pfam accession, forex. (PF02910). Navigate to the domain tab and find the domain architecture that you are interested in and press view scores. That way, only the sequences you are interested in will be downloadable under the tab downloads.
Step-by-Step guide:
1- Go to http://www.ebi.ac.uk/Tools/hmmer/
2- Search with hmmsearch http://www.ebi.ac.uk/Tools/hmmer/search/hmmsearch
3- Use Accession search and type in PF02910
4- Navigate to Domain
5- To the right side of
9742 sequences with domain architecture: FAD_binding_2, Succ_DH_flav_C
Press View scores.
6- It should show (Your results have been filtered) and Query matches (9742)
7- Navigate to Download and download the sequences in the format of interest.
8- Repeat for domain architectures of interest.
I hope you find this helpful.
Kind Regards,
Sara

That does not appear to be possible using PFAM. I would suggest trying a different database. I managed to perform the same query using SMART (http://smart.embl-heidelberg.de/), via entering "Pfam:FAD_binding_2 AND Pfam:Succ_DH_flav_C" under "Domain selection" in the "Architecture Analysis" section of the search screen, which resulted in 9921 protein sequences with those two domains as you can see here.
On the results screen you can select all the sequences and under "Action" at the top of the screen choose "Download Protein Sequences as FASTA file".

It's a good idea but the problem is that it's not specific. I need the proteins that just have these two domains and no more. In these 9921 protein sequences, there are sequences which have other domains that make my search inspecific

you gave me a good idea. although it's inspecific but I could simply make it specific. First in "action selection" part, I chose "download protein sequences as a fasta file" and downloaded the fasta file of all sequences and then I chose "generate a Newick tree and iTol formatted dataset" . Then in "iTOL protein dataset file" part, I clicked on "download as tab delimited plain text" and downloaded all protein identifiers with the name of their pfam domain

Then I used python for deleting protein identifiers which had other domains. After that I could go to the fasta file and delete the fasta sequence of extra proteins


text:
Blueprint replicates
On Encode, it's pretty easy to find out replicate information about the experiments, however I am having some difficulty finding this information on BLUEPRINT DDC, at http://dcc.blueprint-epigenome.eu
Where can I find information which describes technical replicates and biological (isogenic) replicates?
Thanks in advance.

Is that information about the samples not available on GEO?

Blueprint does not have technical replicates, you will find that very few sequencing based experiments ever have technical replicates. As biological replicates: each individual is regarded as a replicate.
This makes sense for the sort of experimental design which blueprint was designed to answer questions such as "Do individuals from group A have higher or lower H3K4me3 than individuals from groups B" or "Does cell type A have higher accessibility than cell type B". In each of these cases the within group variability to be considered is encapsulated in the difference between individuals.
Overall technical variation is assumed to be similar in all samples or at least not correlated with the contrast of interest (this may well not be true in the second example design I gave) it is safe to directly measure it, as long we know its distribution.
As an aside, one of the reasons that ENCODE got away with only two reps is that their reps are effectively technical (repeated preparations from the same cell lines) and thus there is next to no biological variation (which is what our algorithms are set up to find).

Unfortunately this information isn't as clearly laid out in Blueprint as it is in ENCODE. You can view metadata for experiments from one of the many data viewers (IHEC Data Portal, the Blueprint DCC you mentioned, DeepBlue, or the EGA dataset itself), but these tend to not contain technical/biological replicate information like ENCODE does.
You often have to dive into the data from the original publication to get that information


text:
How to get UCSC Genome Browser to show the actual sequence?
Is it possible to make UCSC Genome Browser to show the sequence (nucleotide chain) of the current region?
I fiddled with the various viewing options and couldn't figure it out.

If you zoom in at one specific region, it shows the nucleotide chain on the third line (below the chromosome number). It can't be there by default as the size of default region is more than 100,000bp. You can use the scale on the top of the page to see the size of the region.

In case it’s of any benefit to you, I’ll often check IGV when looking at tracks with that kind of zoomage. Not that it zooms any better than the UCSC browser (although the bases are colour coded!) but you can see your aligned reads individually if you want, along with directional information etc etc. Loads of stuff.

From the top bar:
View -> DNA


text:
Best way to "fish" long reads that match a query sequence
Very simple set-up: We have a PacBio long high-fidelity (HiFi) reads, genome sequencing, and some of those that contain a particular sequence. I want to find out which ones. The length of the query sequence is 7kb, though ideally I would use as query a sub-sequence of only 1.5kb.
The output can the just the read ID or the full read sequence.
We tried DAMAPPER and it gives some decent results but it might not be appropriate to search for the shorter sub-sequence of 1.5kb. It also feels like an overkill for this purpose.
Bear in mind that this is my first foray into long-read sequencing so my mindset is still very much one of short, low error-rate reads.

I have seen the blasr aligner typically used for aligning PacBio reads to a reference.  Is there a reason why blasr would not be a good choice? I guess the very short sequence you are trying to align to might be an issue.
Minimap2 is capable of mapping pacbio reads, and can also do pairwise mapping to see if reads overlap with each other, so this function might help you:
From https://github.com/lh3/minimap2
minimap2 -x ava-pb  reads.fq query_seq.fq > ovlp.paf


text:
How can I colour boxes in Gviz AnnotationTrack in R?
I'm learning the Gviz bioconductor package, I generate a plot as follows:
library(Gviz)
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'))
pdf(file=paste("test","pdf", sep="."))
plotTracks(track)
dev.off()

test.pdf looks like this:
However, what I want is for the boxes to be coloured rather than labelled with the colour word, how to do this?
Note that my real example is much more complex than this.

Did you try the fill argument?
Something like this:
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))

Thanks, seems to work, can't believe I missed that! What colours are accepted as arguments for fill?

Probably all colors, here a list of color names.

@Chris_Rands see ?plotTracks for the options accepted by plotTracks which will then refer you to ?settings for the various things you can tweak (long, long documentation). Finally, availableDisplayPars(track) will show you what you can set for track.


text:
Library for computing BWT-based alignments
I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. As far as I can tell, though, both of these tools and their derivatives are invoked strictly via a command-line interface. Are there any libraries that implement BWT-based read alignment with a C/C++ API?
Python bindings would also be great, but probably too much to expect.

Not C/C++, but just to mention that rust-bio seems to have BWT-related capabilities: https://rust-bio.github.io/

First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.
It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).
Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).
There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:

SeqAn. See the FMIndex tutorial and the Pairwise Sequence Alignment tutorial for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a quick development of a read mapper, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in YARA.
SDSL-Lite. This is a general library for succinct data structures. See the tutorial slides for an idea of how it works. For instance, GramTools are based on SDSL.

Never heard of YARA … better-known mappers also use SeqAn, e.g. both TopHat and Bowtie (though not for the mapping itself, as far as I know).

Sigh, I was hoping to avoid writing BWA's SAM output to disk, but PySAM treats input strings as file names and explicitly doesn't support StringIO objects. Guess I'll have to look into named pipes or something like that.

PySAM supports also file objects – see http://pysam.readthedocs.io/en/latest/api.html#pysam.AlignmentFile. AlignmentFile(filepath_or_object,....: "If filepath_or_object is a python File object, the already opened file will be used."

Yeah, but running proc.communicate gives me BWA's SAM output as a string. It looks like none of the usual string trickery will work, and somehow I'll have to go through the file system to get PySAM to read it.

BWA-MEM can be used as a library. File bwa/example.c shows the basic functionality for single-end mapping. It should give identical mapping to the bwa-mem command line. Header bwa/bwamem.h contains basic documentation. Paired-end mapping is doable, but not well exposed. Several teams, including GATK and MS genomics, have been using bwa-mem as a library.

That is extremely cool.

SeqAn supports BWT tables for use with their parametrisable alignment algorithms.
To use it, follow the general outline for building a SeqAn short-read aligner, and use the FMIndex specialisation instead of — as in the example — the IndexQGram.

If you're looking for a python API, nanopore have created one (although I am not sure how "maintained" it is) https://github.com/nanoporetech/bwapy