text:
Chunk alignment in a name sorted bam for parallel processing
I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.

Could you confirm what you mean by "chunks". Do you mean a genetic locus, e.g. X nucleotides is one "chunk". Alternatively, do you mean subsets of the read depth?

I want to split the alignment into small segments, e.g., 10M alignments per segment. The expected input bam file is name sorted, so roughly 1-10M alignments in the first chunk,10M - 20M alignments in the second chunk and so on. The problem is that I don't want to put the alignments from the same read into different chunks. Thank you.

What kind of parallel processing do you plan? Is is unusual to split by read names instead of loci for parallel processing.

What have you tried? In particular, what python code using pysam have you written to do this? The steps are mostly (A) name sort and (B) iterate over and write a lot of files.

I think it is a bit outside of the scope of this site to write a program for you. Perhaps if you show us what you have tried so far, it would be more possible for us to assist you. In the question's current state, you are asking us to write you a program which incorporates pysam.

Thank you all for your help. Did some search. I think that for fastq files or sam files, it should be easy based on the class here: https://www.blopig.com/blog/2016/08/processing-large-files-using-python-part-duex/
However, for bam files it's more challenging as I don't know how to use 'seek' in pysam. Or maybe just treat bam files as special zipped files and use general gzip reader for this purpose? Or is it possible to do so by using pysam? Thank you all.

Basically, I want to know where all the alignments of a read mapped to, so I decide to process read by read. Here I distinguish reads from alignments.

Which aligner was used? bwa stored alternative mapping positions in the XA tag.

These are RNA-seq data and I use STAR for alignment. I can see four tags, e.g., NH:i:6 HI:i:6 AS:i:49 nM:i:0

You can sort by read name instead of position. This is relevant for whatever you are doing. You should read this thread to understand the bit flags for multiple primary alignments. It should be helpful.

python will be to slow for this job. Here's a awk solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
  BEGIN { part=0; line=n }
  /^@/ {header = header$0"\n"; next;}
  { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" }
    else { print $0 >> part".sam"; line++; }
    last_read = $1;
  }'

The output will be in sam format.
You can controll the number of reads per chunk by setting the n value for the awk command.
EDIT:
A second version that will output in bam:
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
BEGIN { part=0; line=n }
/^@/ {header = header$0"\n"; next;}
{ if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}
  print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam"
  last_read = $1;
  line++;
}

Thank you so much! This works great as well, but I want to use pysam so the code can be integrated into a longer function. I don't even need to save the chunked alignments, but just pass them for parallel processing. It could be quite useful if I want to change my design to split and save a bam file into smaller bam files. I think that many other people may find this solution helpful as well.

I was considering awk too. The problem with Perl/Python is the code will process the sequence data "[horizontal] line by [horizontal] line", when awk can do this "[vertical] line by [vertical] line "

I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing code if going the be the slow bit.
import pysam
infile = pysam.AlignmentFile("input.bam")

chunk_size = 10000000
outfile_pattern = "output_segement%d.bam"

chunk = 0
reads_in_this_chunk = 0
old_name = None
outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

for read in infile.fetch(until_eof=True):

    if old_name != read.query_name and reads_in_this_chunk > chunk_size:
        reads_in_this_chunk = 0
        chunk += 1
        outfile.close()
        outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

    outfile.write(read)
    old_name = read.query_name
    reads_in_this_chunk += 1

outfile.close()

Thank you, this works great! What I can think of to optimize the code a little bit is 1) I don't need to check the name of each alignment. Only after reading a specified number of alignments, I check the read name, read again until I find a read with a different name. 2) It seems un-doable for pysam to pass the start and stop positions of chunks to multiprocessing.pool.

Not sure how can I accept this answer.

@jenny00 the comparison of read names will take next to no time. The slow bit is reading the read into python, once its in manipulations are fast.
And as for 2: No, pysam is not compatible with python's simple multiprocessing API because you can only pass picklable objects through to the child processes and AlignedSegments (or reads) are binary objects - they are actually a thin adaptor layer on the C++ based htslib object.

Thanks, I just noticed that. So I have to put the necessary information from each alignment into a dictionary and pass dictionaries down for parallel processing. It does not save much time compared to just using one process.

If you are using pythons multiprocessing.Pool infrastructure, you might find that the overhead associated is high - each submission to a Pool causes the creation of a new python process, which does its processing and then the process is closed down (i think). You might like to consider the more complex to implement asyncronous worker and processing queue type set up. Also try increasing the size of the chunks - the fewer chunks, the small the number of process spin ups and close downs.

Thanks. Increasing the chunk size and using map_async help. But it is only twice faster using four cores, compared to using a single core.


text:
Finding the members in a confusion matrix
I have inferred a confusion matrix of training and test set by neural network. I want to know which members are in the confusion matrix.
> table(test$TRG,predicted.nn.values$net.result)

    0 1
  0 5 4
  1 5 3
>

> head(train[,1:4])
  TRG      ALB     AQP9   CALML5
1   1 5.865827 8.190945 6.705303
2   1 6.998435 8.424261 8.505591
3   1 7.424512 8.716471 7.249556
4   1 7.442049 8.325263 8.809286
6   1 5.893411 8.199990 6.677618
7   1 7.288030 9.143510 7.088598
>

> head(test[,1:4])
   TRG      ALB     AQP9    CALML5
5    1 6.369307 7.954310  6.920290
8    1 6.181902 8.651442  7.225389
10   1 6.119359 9.345270  6.829623
13   1 7.775533 9.016272  7.976813
14   1 5.913656 9.484457 10.609013
18   1 6.603138 7.908560  8.827173
>

How can I know which patients are false positive or true negative? For example I have 5 patients predicted correctly, how can I know the name of these patients?

Could you please clarify what is each object and how was it obtained? Without knowing the details of the object we can only guess. What do you mean by "which members are in the confusion matrix"?

This does not look a great confusion matrix, but it is a very, very cool approach. The cool thing about machine learning is the ability to examine lots of models within confusion matrices.

To answer your question, you want a heat map of the regression weights to establish who is who, i.e. which patient is false positive etc... This is the only way I'm aware of to analyse a neural network. This looks like R, I only know Scikit learn.

Sorry, why you are mentioning my approach is a very, very cool whereas I better change the model? Which aspect is cool then?

I have changed my answer to be more diplomatic, apologies. Neural networks are very cool per se largely because of Google - albeit they are champions of reinforced learning. The problem with neural nets is the difficulty in accurately understanding what each component (patient in this case) has contributed to the overall result and why. Other approaches are much clearer on this. Ultimately they are just one approach amongst many machine learning approaches. I like the whole area - except for having to learn Python.

source: Stack Exchange

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

sort: votes


text:
Which gene I should select from this qqplot
I have a qqplot of my whole genome sequencing data; A plot is for showing possibly significant driver genes. I tried to read about qqplot though but people only say about the skewedness while I want to know from these two genes which one are more likely to be driver

Given that TP53 is the most significant and is already known to have driver mutations in cancer it would seem to be the logical choice.

Thank you, does that mean I must ignore MAP1S as a candidate even though it has significant p-value but non-significant q-value?

Nothing I wrote suggests that you must ignore MAP1S, your question indicated you only wanted a single candidate. If you're fine with more then take more.

You can't know which is the driver and which is the carrier. At most you can say that a specific gene deviate more of the expected underlying hypothesis. See also other resources online.
You also seem to ignore other genes that deviate more of your null hypothesis. I recommend to plot the histogram of the p-values to see if the distribution is uniform

I think it's pretty clear from the QQ-plot that these p-values are uniformly distributed (except for the twenty or so in the tail of the distribution)?

@winni2k Usually the QQ-plot is centered around 0 (if using qqnorm), in this plot it is not centred, also a expected and observed pvalue above 1 indicates that this is not a normal QQ-plot. That's why I recommended the histogram.

Ah! I see. I think the x and y-axes must be -log10(p-value) as shown for example here, and not quantiles as would be reported by qqnorm. #missingaxislabels

This looks to be a pretty good Q-Q plot. My recollection of doing Q-Q plots is that a good Q-Q requires a "more or less" linear relationship for the model to be considered okay. Again this model looks good, its 1:1, albeit it does deviate a little (for the gene of interest). TP15 falls outside the 1:1.
@Devon Ryan states TP15 has good biological credentials, which is clearly important. It maybe worth checking  MAP15 because this is "more or less" within your model (for my Q-Q plots definitely within the model :-) ).
In summary, its just my opinion - that on the basis of the Q-Q plot the driver is the not the expected candidate.

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


text:
How do I perform a pathway related grouping of genes?
I have a list of gene names from the analysis of my mouse single seq data. I want to perform a simple grouping of the genes based on the pathways to which they belong or pathways through which they relate. Results from my online searches strongly recommend the ingenuity pathway analysis (IPA) tool. My laboratory has a license for the tool, but I need help with the workflow/steps. I have no prior experience with pathway analysis, and I have a background in R and Python programming. Can anyone also suggest any link which can help me proceed with the required IPA analysis? I will also be very grateful for any other tools (databases, R or Python packages) that can be helpful here.

Have you tried watching any of the many videos Qiagen has made for how to use IPA? They're quite helpful. Note that you'll want to upload a list of gene identifiers with fold-changes (and possibly p-values, unless you've already subset them).

What are you interested in? The title ask one question but the body ask another one. Which is your biological question of interest? Could you also explain what steps have you done to find help with the IPA software?

FYI: Question is also posted on biostars.

The DAVID tool is always interesting, and easy to try: https://david.ncifcrf.gov/gene2gene.jsp
You submit a list of genes, it uses several different annotation databases to annotate the gene list, and then it clusters the results and tests for statistical enrichment. You can potentially learn about functional enrichment in your gene list from several different independent annotations.


text:
Determining Read Groups
Which Read Groups are correct:
java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84

or
java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84_S23_L005 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84_S23_L005

Or is a better way to determine read groups?
Thank you in advance,

Usually it doesn't actually matter, but if you want it to be very correct:
RGID=SomeLibraryID \
RGLB=SomeLibraryID \
RGPL=illumina \
RGPU=Some_Sequencer_ID \
RGSM=SomeSampleID

That assumes you have one library per sample and aren't splitting the output BAM files by lane.

Is it possible to get it automatically determine by reading the FASTQ header or filename?

If you have control of how the data is produced then yes, you can set all of that automatically. Otherwise you can't be certain the file/read names are meaningful.


text:
Alignment with arbitrary number of mismatches or gaps
I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (<7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:
echo "GGAGAACAAGTGGCAGATAGAGG" | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4

I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:
0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU

I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.
The last three bp will have to match exactly - I was going to enforce this manually later on, maybe there's a better option for doing this. Taking the indels into account is also a secondary aim.
I'm fine with ~10k alignments per query for <=6 mismatches plus the last 3bp matching exactly. Haven't done the combinatorics but that sounds realistic to me - question is how to make bowtie do it, or a different tool.
What would be the right way of achieving this using bowtie, or a different tool?

This is not something mappers will be good at, because you're going to hit most of the genome. Bowtie won't do it, and I don't think LAST will either. The majority of fast mapping programs rely on seed-based mapping where the seed is essentially free of errors.
Doing only mismatches makes it much easier to determine, but it's still too relaxed. Consider that the example sequence you have provided is 10 mismatches away from GGGGGGGGGGGGGGGGAGAGAGG, and also 10 mismatches away from AAAAAAAAAAAAAAAGAGAGAGG (with the last 3 bp matching exactly)
Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:
for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp

Seems like I don't have enough reputation to upvote, but thanks!

I've ended up using a package called BatMis for this, which can detect alignments with an arbitrary number of mismatches (not sure about indels). It is quite performant and yields about the order of magnitude of alignments that I was expecting (several thousands for alignment of a 20bp read to the whole human genome with up to 6 mismatches).

Did you see

-k  Report up to  valid alignments per read or pair (default: 1).

Did see that, but I thought --all (or -a in short) means it outputs all found alignments, ie supersedes -k? Using -k 5 instead of --all still only yields one result.


text:
how to obtain only one cluster's data (.csv file) in Seurat
I was using Seurat to analysis single-cell RNA Seq. And I was interested in only one cluster by using the Seurat. Does anyone know how to achieve the cluster's data(.csv file) by using Seurat or any other method that can do this function? Thank you!

A cluster may contain different types of data. Do you want the specific cells, expression for each cell, average expression, cluster markers, etc.?

Yes even a cluster contain different types of data. I still  want to obtain that cluster's data including its each cell and gene expression for further analysis.

I think that the Seurat tutorials had some information about this. Did you read them? Have you tried something?

You can access data within the Seurat object using GetAssayData, and extract a list of cell names for the cluster you're interested in using WhichCells:
library(Seurat)

cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)

Hi TimStuart, thank you for your advice, and it works well in the RStudio. But when I use this script in cypress it shows the error "Error in as(Class = "matrix", object = expr), could not find function "as"
Execution halted". Do you know how to solve this problem?

I don't know what cypress is, but if you're running it in a script you'll need to load the methods and utils libraries

Hi TimStuart, thank you for your reply, and if it is possible could you please tell my what kind of the utils libraries that I should load?

library(methods); library(utils)

Hi TimStuart, nowadays thank you for the Seurat's update, we can analysis several samples together. And I would like to know whether we can obtain the specific sample's specific cluster gene expression profile by using the same code "cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)"?


text:
Strategies to disentangle the effect of two genes on phenotype
I am analyzing the effect of one putative tumor suppressor (and its loss) on a specific cancer phenotype. This putative tumor suppressor is known in literature to have an effect on the growth of tumors in mouse models.
However, in the same region of this gene there's a second well-known tumor suppressor - actually a very powerful one. By looking at the data from cBio portal, there's a very strong co-occurrence of loss on both tumor suppressors in a pan-cancer analysis, that is they are almost always lost together.
Therefore, I would like to disentangle the effect of the two genes, and to analyze only the effect of my putative candidate on the phenotype. The simplest approach is to use only the samples that have a single loss (only the loss of my putative candidate, but not the loss of the confounding tumor suppressor). However, only ~10% of all samples shows such behavior, and this 10% is distributed across all cancer types, decreasing the statistical power of the analysis - since I would like to test each cancer type separately to avoid mixing incompatible data.
Which kind of strategy/analysis can be used in this case?

How do you model the effect of the genes? I would use all the data (or cancer specific data) and indicate which sample have lost the gene of interest, the other or both genes. That way you could disentangle the effect of each gene and the interaction of both

Maybe you could fit a mixed model to predict phenotype, using the two genes as fixed effects and cancer type as a random effect, and ask the model for the separate effects of each gene, and the interaction between them.

Or even more simply using generalized linear modeling (GLM), where the regression weights are the genes and the accuracy of the model is measured in the residual. You also Q-Q plots to assess the fit of the model to the data.
If the residual is high and Q-Q "wobbly" transformations works wonders.

Machine learning such as a neural network is a very cool way to do this calculation, but IMO you need Python - although people on this board do this in R. The concept is very easy you are catagorising (keep in mind I'm not an oncologist) gene vs. cancer resistance and you want neural network, linear-based regression etc.. to make the call (you will see it visually). You then use the weights to assess the relative contribution of each gene.
The 10% issue is formally assessed, its not hard to obtain the power because its a common problem. The formal way to assess this, from memory, is you assign a "null model" and the accuracy of a ML model needs to a lot higher than the null model. I forget the details, I've never been in these situations. There's a specific technical name (which I should look up).

The explicit GLM / GLMM (non-ML) is the easiest approach however. If you use a Mac there is a great (relatively cheap) package for GLM, but any good software package will do. I assume GLMM is available in R (certainly Stata).


text:
Help interpreting my phylogenetic tree construction of bacterial species in same genus
I constructed a phylogenetic tree using mole blast, and it created a neighbor joining tree, so I used the alignment of those sequences with mine (1-43) in mega with TN93+G+I model but I don't understand how to interpret why my species in question (29) is so much different (or maybe novel?).
Attached is my original tree (Picture 1) and my bootstrap analysis (Picture 2):

Are these the same tree with different roots? The first tree is probably rooted by middle point. How did you root the second tree?

Yes that is correct. The interpretation is easy and relates to the population structure of two definite bacterial species within the same genus. I will provide a detailed interpretation shortly.

The top phylogeny is likely mid-point rooted, but I suspect the rooting is fine. The lower diagram has been rooted on the query sequence. It is taxonomic analysis using 16S which in bacteria is considered to represent species delineation at a given percentage of divergence. I thought it was 5%, which is impossible based on the phylogram in the top picture.
Bootstrap values lower than 80% mean that the clade is not supported (monophyletic group is a better term). The delineation between 10001 and 10013 against the massive ingroup of 13 16S sequences is robustly supported by bootstraps of 95% and 99%. The divergence between them is around 2%. There is no taxonomic distinction between the 13 16S sequences based on the phylogram.
One question is whether this is full length 16S? If it is there appears to be a species designation confusion in the lower clade supporting 10001 and 10013 because it is identical to to Dermococcus species and there is no way to assess which one is valid, nor whether they are really members of the same species. The very small divergence suggesting in D. abyssis might also call this into question. 10013 is in fact a robust outgroup to 10001 and the three Dermococcus spp. but the genetic distance involved is small.
However, this is not the major result. You have recovered what appears to be a large population (if that is how you sampled it) with close genetic homology to Dermococcus strain DSM20448. Given the complexity of species lineation within this genus it is impossible to say whether this represents the same species as DSM20448 but you could call it a "species complex".
Overall the population structure (if that is how it was sampled) is heavily skewed towards a bacterial species which is hitherto received minimal attention. This population skew (if that is how it was sampled), could be reversed in different ecotopes/niches/ecoregions/patients?
In summary, I wouldn't get caught up in systematics because in this genus its a bit messy. I would focus on the population structure and these results are very clear.

Thank you so much! It is only a partial sequence so the bootstrap being low isn't a surprise.


text:
Combat for multiplatform batch correction
I have been looking ways to use data from various micro-array platform such as agilient, rosetta/merk. Some of the established method I came across is combat, but I'm not sure if it can be done for cross platform.
A similar question here where the answer is kind of can't be or shouldn't be done.
Any suggestion how to take account of batch effects across multiple platforms.

Usually with microarrays you want to make a case/control comparison, so I am going to assume that.
Data from different array platforms is generally difficult to compare:  each platform is measuring potentially a different part of the expression of the gene (different exons or different regions), each platform is likely to need different normalisation approaches, and have different dynamic range.
If you have case/control data from several different platforms, then the most robust thing to do is to analyse each platform separately, and produce a list of up and down-regulated genes according to each platform.
Then, it would be reasonably safe to combine the lists of genes in a 'meta-analysis' approach.  I would expect that in a lot of cases, 'up' and 'down' regulation would be correlated across platforms, but I would not necessarily expect very strong correlation in the magnitude of changes between platforms for individual genes.

so are you saying its better not to use combat or any other algorithm for that purpose? due to multi platform issue , we have out of one control which for other comparison so meta analysis may not not be possible to different patform

You could try combat, I don't have experience with that.  In general I would say it is difficult/not advised to mix platforms in a single analysis.

okay .will wait for other answers if anyone has done something similar..to get their views..


text:
GEOquery errors due to list input
I want to extract the expression matrix from the GSE object in GEOquery.
library(GEOquery)

gse1 <- getGEO("GSE97168",GSEMatrix=TRUE)
datExpr = exprs (gse1[[1]])
datExpr

GSM2553017 GSM2553018 GSM2553019 GSM2553020 GSM2553021 GSM2553022 GSM2553023 GSM2553024 GSM2553025
     GSM2553026

Table(GSMList(gse1)[[1]])[1:5,]

Table(GSMList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GSMList’ for signature ‘"list"’

Table(GPLList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GPLList’ for signature ‘"list"’

Meta(gse1)$description
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Meta’ for signature ‘"list"’

Columns(gse1)
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Columns’ for signature ‘"list"’

Why it's not working with R version 3.5.1?

Welcome to the site. Are you sure this should be working? In the vignette Table is only used to "Look at data associated with the GSM" not to retrieve the data expression data or information about a specific GSM

Use GSEMatrix=FALSE, the information you want appears to be in the soft file, rather than the matrix. Yes, this is moderately annoying since it means you need two objects to describe this particular dataset.


text:
Where to Download Cancer Raw Reads (fastq)?
Does anyone know where I can download cancer raw reads (fastq files), tumor and Germline for non-humans?
I wanted to make a study with human data but I don't have the control access to download raw reads.

So do you want human or non-human? Your question is confusing. For human, have a look at Is there a publicly available tumor-normal sample?. Is that what you need?

Please tell us when you cross-post your question. This question is also asked on biostars.

You can try to browse SRA (sequence read archive) at https://www.ncbi.nlm.nih.gov/sra.
To download fastq files from SRA, you'll need their sratoolkit: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software.
The easy command to download fastq files by their accession number is :
fastq-dump --split-files accession_number


text:
Error in adding data in columns of R data frame
This is my file
df1 = read.csv("HSC_LSC_BLAST_karyotyope.txt",header = TRUE,sep = ",",row.names = 1)

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- rep(NA, ncol(df1))
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
#metadata$Group[seq(13,16,1)] <- 'Mono'
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'

>s metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy' Error in `$<-.data.frame`(`*tmp*`, Cytogenetics, value = c("Healthy",  :
> replacement has 4 rows, data has 23

Im trying to add other metadata information based on clinical information apart from giving sample labeling , im certainly doing something wrong but im not sure what it is. As my idea is the another column named Cytogenetics which would have respective information status but im getting that error.

I know this is not really an answer, but I strongly discourage you from doing what you do, it's dangerous to impute data manually because it's prone to human mistakes - that you might mess up which columns should carry with metadata.
Actually, I think I have spotted a mistake already. This is your table using the corrected version of your code (with metadata$Cytogenetics <- NA, see Devon's answer):
>metadata
        Group Cytogenetics
HSC1      HSC      Healthy
HSC2      HSC      Healthy
HSC3      HSC      Healthy
HSC4      HSC      Healthy
Blast11 Blast           NK
Blast12 Blast           NK
Blast1  Blast           NK
Blast2  Blast           NK
Blast3  Blast           NK
Blast4  Blast           NK
Blast6  Blast           NK
Blast7  Blast           NK
Blast8  Blast           NK
Blast9  Blast           NK
LSC1    Blast           NK
LSC2      LSC           NK
LSC3      LSC           NK
LSC4      LSC           NK
LSC6      LSC           NK
LSC7      LSC           NK
LSC8      LSC           NK
Blast5    LSC     Abnormal
LSC5      LSC     Abnormal

I sort of think that the column LSC1 should be in the LSC group and visa versa the Blast5 should be in Blast group. In this case it seems that it's derived from name, you can actually write a code that will extract the group from the name.
df1 = read.csv("HSC_LSC_BLAST_karyotyope.txt", header = TRUE, sep = ",", row.names = 1)

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- sub("[^[:alpha:]]+", "", (colnames(df1)))

However, I don't know how do you know which sample is Healthy, NK or Abnormal. I suppose you have a table where this is written, if it's the case, use that table, don't impute the data manually.

yes actually the data im taking is realy big so im going for selected sample as you said to do it manually its prone to error ,one more question i want to add respective blood profiling data which are counts such as WBC and other counts, do i have to just make another data frame and add ,or i have to make a column of counts and add it as you shown?

"I don't know how do you know which sample is Healthy, NK or Abnormal. I suppose you have a table where this is written" yes they have categorized based on karyotyping as well as other information

This is not really the good platform for discussing good practices. I would just load the two tables and use function merge to aggregate the data together.

Add the Cytogenetics columns and then assign things to a subset of it.
metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- NA     # you don't need to write a vector; just NA does the job
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
metadata$Cytogenetics <- NA    # initiate the column with NA values
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'


text:
Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3
I downloaded the annotation of the C. elegans genome in GFF3 format from Ensembl.
I typed the following command, hoping to get the header of the file (lines starting with #).
grep '^#' Caenorhabditis_elegans.WBcel235.95.gff3

This returned 46798 lines, many many more than the dozen or so header lines. After searching the file for the # character, I noticed many lines that contain only three hashtag (#) symbols.

Why are these lines there?

Could you provide the lines and the context around the lines?

I added a screenshot. There are too many lines to show them all here. There is no constant interval between the ### lines.

Next time, please don't post images of text. Just copy the text directly into your question and use the formatting tools. images are harder to search, make the page heavier to load (important for users with slow connection) and cannot be copied (I can't extract your file to demonstrate a command that solves your problem).

There are some good answers so far, but I don't think any of them fully communicate the significance of the ### directive. The GFF3 specification states:

This directive (three # signs in a row) indicates that all forward references to feature IDs that have been seen to this point have been resolved. After seeing this directive, a program that is processing the file serially can close off any open objects that it has created and return them, thereby allowing iterative access to the file. Otherwise, software cannot know that a feature has been fully populated by its subfeatures until the end of the file has been reached. It is recommended that complex features, such as the canonical gene, be terminated with the ### notation.

In other words, without the ### directive you have to load the entire GFF3 file into memory to make sure you've inferred all feature relationships completely (as encoded by ID and Parent attributes). With the ### directive, it's possible to process the annotations in smaller chunks, with a MUCH smaller peak memory usage.

So, you don't want to get rid of them and process the file as a single object.

@MichaelG Right-o! Note, that having ### directives only makes it possible to process GFF3 files efficiently. It is not a guarantee that any particular tool will utilize these directives.

The three # are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the #.
To grep just the header, which has 1 or 2 # from the line start on, you can use extended regular expression:
$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

Which means: "grep all lines that started with 1 or 2 # and are not followed by #"
If you like to create a new gff3 to get rid of this separator lines you can do this:
$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3

The -v prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")

I know "negator" is made up ... The technical term is "negation operator", (https://en.wikipedia.org/wiki/Negation). awk is more powerful than sed, as finswimmer demonstrates. ... original post edited

I discussed with  @terdon. I didn't realise sed/grep had negation operator [^#], I thought this was regex such as Perl..  ... so .... grep "^#.[^#]" ./Caenorhabditis_elegans.WBcel235.95.gff3  ... is another solution

As others have stated, those are just there to separate the entries for easier parsing. They enable you to do nifty tricks like:
$ awk -v RF='###' '/Y74C9A.3/' Caenorhabditis_elegans.WBcel235.95.gff3
I   WormBase    mRNA    4116    10230   .   -   .   ID=transcript:Y74C9A.3;Parent=gene:WBGene00022277;Name=Y74C9A.3;biotype=protein_coding;transcript_id=Y74C9A.3
I   WormBase    three_prime_UTR 4116    4220    .   -   .   Parent=transcript:Y74C9A.3
I   WormBase    exon    4116    4358    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e5;constitutive=1;ensembl_end_phase=-1;ensembl_phase=0;exon_id=Y74C9A.3.e5;rank=5;version=1
I   WormBase    CDS 4221    4358    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    5195    5296    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e4;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e4;rank=4;version=1
I   WormBase    CDS 5195    5296    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    6037    6327    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e3;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e3;rank=3;version=1
I   WormBase    CDS 6037    6327    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    9727    9846    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e2;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e2;rank=2;version=1
I   WormBase    CDS 9727    9846    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    CDS 10095   10148   .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    10095   10230   .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e1;constitutive=1;ensembl_end_phase=0;ensembl_phase=-1;exon_id=Y74C9A.3.e1;rank=1;version=1
I   WormBase    five_prime_UTR  10149   10230   .   -   .   Parent=transcript:Y74C9A.3

To extract the headers, you just need the lines that start with exactly 2 #:
$ grep '^##[^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942

Or, to also include the comments, lines that start with a #, then either a # or a ! and then a non-#:
$ grep '^#[#!][^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06

Or just a #, any character and then a non-#:
$ grep '^#.[^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06

@charlesdarwin it prints the output I show in my answer :) Basically, it sets the string ### as the record ("line") separator and then print the record that matches Y74C9A.3. The default action when something evaluates to true in awk is to print the line. So, /Y74C9A.3/ will be true on lines that contain Y74C9A.3, and since we have used ### as the line separator (record separator) it will print the entire entry for Y74C9A.3.

@charlesdarwin why did you add the space after the -v? Do you know of an awk implementation that requires it? I fear that makes it look like a "verbose" flag.

my awk complained without the space. awk: invalid -v option. I am on macOS X high sierra. using awk version 20070501. crazy number there! is it a date?

@charlesdarwin ah, yes, I tested on GNU awk. Thanks! And yes, that looks like the 1st of May, 2007.

Hmmm, I dunno. If I was doing this I'd want an easy way to parse the file and ### would be a good marker. Anyway just pop a space character after the #. Using sed ...
sed -n '/^# /p' < ./Caenorhabditis_elegans.WBcel235.95.gff3

That should remove the ### lines
Also head -20 and sed -n '/1,20/p' ./Caenorhabditis_elegans.WBcel235.95.gff3 would get around the problem with a bit of tweaking.

The grep verion is (removing unnecessary jargon)
grep '^# ' Caenorhabditis_elegans.WBcel235.95.gff3

This is the same command,
grep '^#\s' Caenorhabditis_elegans.WBcel235.95.gff3

The point behind the code is there is no negation operator, in sed (I assume grep is the same). In Perl you would simple [^#]{3} (need checking) in your regex (remove all ###). In sed (or grep) you can exclude them using the ^ (beginning of a line) and suffixed with a "space" with only a single character in the middle all 3 letter patterns are removed.

Another approach without the space character is
sed -n '/^###/d; /^#/p' ./Caenorhabditis_elegans.WBcel235.95.gff3

Albeit this may require,
sed -i '/^###/d; /^#/p' ./file

what is your sed command doing? There is no output and no change to the file when I run it.

I've edited the post to include the command for your grep approach and included an explanation. It is worth learning sed BTW, it has some nice tricks, but it is not as important as awk.

Your sed command does nothing. I think you meant to use -i there, but I don't really see how. There are no lines matching ^#  in the file. To remove the ### lines, you would want something like sed -i '/^###/d' file.gff. And you can negate character classes in both sed and grep just as you do in perl: [^#] means "not #". This is a general regex feature available in all common regex flavors. It isn't a PCRE thing (although GNU grep can deal with PCRE as well, just use the -P flag). But note that [^###] does not mean "not ###", it just means "not #".

sed -i '/^###/d' would not work BTW, becaue sed -n '/^###/d' does not. sed '/^###/d' is an improvement. One way to get it working would be >sed -n '/^##/d; /^#/p' ./file ... I've edited the post accordingly

I suspect sed -n has been deprecated with implementations of sed -i ... sed -n is discussed here https://stackoverflow.com/questions/6255796/how-the-n-command-works-in-sed … ...... It appears sed -i is not supported in all sed implementations, https://stackoverflow.com/questions/18527365/what-does-sed-i-option-do .... SED (-n): is direct from O'Reilly's famous "sed & awk" Doughty & Robbins, e.g p71 ... sed -n works on OS X and Susie [bash and tcsh]. PERL: [^###] is identical to [^#] but [^#]{3} I think would resolve the issue.

No, sed -n is in no way deprecated, it's defined by POSIX.  And it has nothing to do with -i ("make the changes in the original file"), -n means "don't print unless I tell you to". And -i it is supported (with minor format differences) in all sed implementations you are likely to use as a bionformatician (GNU sed on Linux, BSD sed on BSD and OSX, even busybox sed on embedded systems). Anyway, sed -i '/^###/d'  works just fine, try it out! It's just that it edits the original file, so you don't see any output.

By the way,  sed -n '/^###/d' also "works", it just doesn't do anything useful. That means "delete lines starting with ### but in any case, don't ever print anything". Perhaps you meant something like sed -n '/^###/!p'.

There must by different versions of sed. I do not have any negation operators on mine - this was always my understanding of its weakness, ie your '/^###/!p'. I do agree your /^###/d' can be stacked into a /^#/p command and act as a work-around for negation, which I'd never thought about.

Let us continue this discussion in chat.


text:
Making a bed file for RSeQC
I making a bed file for RSeQC, so it can do things like compute the number of reads from exons, introns, 5"UTRs, etc.
I want to use a bed file that corresponds to my GTF file, so I use gtf2bed to make a bed file, like this:
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' Homo_sapiens.GRCh38.gtf | gtf2bed - > $Homo_sapiens.GRCh38.bed

This file throws some errors, among other when using it with junction_annotation.py:
Reading reference bed file:  /rst1/2017-0205_illuminaseq/data/references/Reference_Genomes/GRCh38.95/Homo_sapiens.GRCh38.87.bed  ... Traceback (most recent call last):
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 125, in <module>
    main()
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/bin/junction_annotation.py", line 109, in main
    obj.annotate_junction(outfile=options.output_prefix,refgene=options.ref_gene_model,min_intron=options.min_intron, q_cut = options.map_qual)
  File "/rst1/2017-0205_illuminaseq/tools/snakemake_conda_env/f4c16777/lib/python2.7/site-packages/qcmodule/SAM.py", line 3762, in annotate_junction
    exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
ValueError: invalid literal for int() with base 10: 'gene_version'

and when I compare this bed file to the RSeQC provide (dated) one these are the differences:
My bed file I created:
1       11868   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00002234944"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       11868   14409   ENSG00000223972 .       +       havana  gene    .       gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_id "";
1       11868   14409   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; tag "basic"; transcript_support_level "1";
1       12009   12057   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "1"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001948541"; exon_version "1"; tag "basic"; transcript_support_level "NA";
1       12009   13670   ENSG00000223972 .       +       havana  transcript      .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; tag "basic"; transcript_support_level "NA";
1       12178   12227   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001671638"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12697   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "3"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001758273"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       12612   12721   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "processed_transcript"; exon_id "ENSE00003582793"; exon_version "1"; tag "basic"; transcript_support_level "1";
1       12974   13052   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "4"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001799933"; exon_version "2"; tag "basic"; transcript_support_level "NA";
1       13220   13374   ENSG00000223972 .       +       havana  exon    .       gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000450305"; transcript_version "2"; exon_number "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-201"; transcript_source "havana"; transcript_biotype "transcribed_unprocessed_pseudogene"; exon_id "ENSE00001746346"; exon_version "2"; tag "basic"; transcript_support_level "NA";

The RSeQC provided bed file:
1       11868   14409   ENST00000456328 0       +       11868   14409   0       3       359,109,1189,   0,744,1352,     DDX11L1
1       12009   13670   ENST00000450305 0       +       12009   13670   0       6       48,49,85,78,154,218,    0,169,603,965,1211,1443,        DDX11L1
1       17368   17436   ENST00000619216 0       -       17368   17436   0       1       68,     0,      MIR6859-1
1       14403   29570   ENST00000488147 0       -       14403   29570   0       11      98,34,152,159,198,136,137,147,99,154,37,        0,601,1392,2203,2454,2829,3202,3511,3864,10334,15130,   WASH7P
1       29553   31097   ENST00000473358 0       +       29553   31097   0       3       486,104,122,    0,1010,1422,    MIR1302-2
1       30266   31109   ENST00000469289 0       +       30266   31109   0       2       401,134,        0,709,  MIR1302-2
1       30365   30503   ENST00000607096 0       +       30365   30503   0       1       138,    0,      MIR1302-2
1       34553   36081   ENST00000417324 0       -       34553   36081   0       3       621,205,361,    0,723,1167,     FAM138A
1       35244   36073   ENST00000461467 0       -       35244   36073   0       2       237,353,        0,476,  FAM138A
1       52472   53312   ENST00000606857 0       +       52472   53312   0       1       840,    0,      OR4G4P

Why are the files so different? How would I make a simpler bed file form my GTF file?

Please show us an excerpt of your GTF file (a few lines) and i) the output you get from it and ii) the output you want to get from it. As it stands, I don't understand what you mean by "simpler bed file" nor why the files are different (can't know that if we don't know what data you are starting from).

UCSC Utilities
Such output results when using the UCSC utilities gtfToGenePred and genePredToBed in series. Builds for macOS (x86_64) and Linux are available directly from UCSC. Alternatively, Conda builds are available through the Bioconda channel, under the package names ucsc-gtftogenepred and ucsc-genepredtobed.
Example: Homo sapiens (Ensembl Release 106)
## (optional install through Conda)
#mamba install -c conda-forge -c bioconda ucsc-gtftogenepred ucsc-genepredtobed

## download GTF
wget ftp://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.gtf.gz

gzip -cd Homo_sapiens.GRCh38.106.gtf.gz |\
  gtfToGenePred /dev/stdin /dev/stdout |\
  genePredToBed /dev/stdin /dev/stdout |\
  head

Output BED
1   1211339 1214153 ENST00000379236 0   -   1211554 1214127 0   7   286,129,197,67,102,123,171, 0,364,602,1298,1652,2323,2643,
1   1211339 1214138 ENST00000497869 0   -   1214138 1214138 0   5   493,197,67,794,156, 0,602,1298,1652,2643,
1   1212018 1213498 ENST00000453580 0   -   1213498 1213498 0   4   120,67,102,104, 0,619,973,1376,
1   1203507 1206571 ENST00000328596 0   -   1203590 1206571 0   4   453,88,123,187, 0,891,1862,2877,
1   1203507 1206592 ENST00000379268 0   -   1203843 1206571 0   5   461,203,88,123,208, 0,526,891,1862,2877,
1   1203843 1205680 ENST00000486728 0   -   1203843 1205463 0   4   125,203,88,311, 0,190,555,1526,
1   1203843 1206571 ENST00000379265 0   -   1203843 1206571 0   5   125,182,88,123,187, 0,211,555,1526,2541,
1   1471764 1497848 ENST00000673477 0   +   1471884 1495817 0   16  325,77,102,60,70,166,70,156,57,126,125,52,71,168,109,2364,  0,5509,6879,7284,9102,10373,10780,13251,14017,14345,14779,16098,17439,18492,18798,23720,
1   1478025 1497848 ENST00000472194 0   +   1497848 1497848 0   14  720,60,70,166,70,156,57,126,125,52,71,168,109,2364, 0,1023,2841,4112,4519,6990,7756,8084,8518,9837,11178,12231,12537,17459,
1   1479048 1482662 ENST00000378736 0   +   1482662 1482662 0   4   60,70,166,118,  0,1818,3089,3496,

Thanx for this answer, I'm not really in a position to check this quickly anymore, but I'll label this as solved.


text:
BAM to gene expression matrix (UMI counts per gene per cell),10X
I am trying to reproduce some results of a scRNASeq experiment. However I am new to the server-side aspect of such analyses and am very confused at the moment.
The data provided by the authors of the paper is in .BAM format and from there I wish to derive a gene expression matrix (UMI counts per gene per cell). The authors stated that they did this using the 10X genomics pipeline, however it seems that .BAM outputs are sort of the "final version" of the pipeline and, the pipeline only really deals with truly raw sequencing files (bcl/fastq). I am trying to learn how to do all this processing myself as I will need to eventually. Unfortunately, most of the documents I have read end up making me more confused. I know BAM files are technically already processed and aligned but, how would I get the counts from the file?
I have considered converting it back to fastq format and following the pipeline from the beginning, but it just seems like I would be moving backwards by doing such.
I have also found some R packages that can integrate BAM files in the environment like Rsamtools, but I would much rather do this on the server and then upload a csv into R for the downstream analysis as the files are +60GB each.
Essentially I am just asking for some advice, or a vignette/article I can read that would clear up how to proceed with the BAM outputs to derive the gene expression data. I have already read the cellranger pages on the 10X genomics website but they were of no help in answering this quesion.
Thank you for your time,

How to re-analyze 10X BAM files?
This is a great question and honestly, I don't think there was an easy way to do this at the time the question was asked. The reason is that if you want to re-do the authors' analysis to get the gene expression (gene-cell / feature-barcode) matrix and they used cellranger (the official 10X pipeline software) to do that, you will need FASTQ files as input - because "10x pipelines require sequencer FASTQs (with embedded barcodes) as input."
Therefore, to get FASTQ files from 10X BAM files, there's a relatively new tool provided by 10X Genomics: bamtofastq. Then you can simply use those FASTQ files as input to cellranger.
More info here: https://support.10xgenomics.com/docs/bamtofastq
How to correctly download 10X BAM files from NCBI?
The only caveat is that you need a BAM file generated directly by 10X's cellranger (or the respective 10X pipeline, if not dealing with gene expression) - that means that a BAM file obtained by downloading an SRA from NCBI and converting to BAM won't work; you need to get the original BAM file directly (often found among the originally submitted files, under "Original format").
To give you some additional help here (because figuring this out is also not trivial): an original 10X BAM file (if originally submitted in this format) can be downloaded from NCBI using prefetch from sra-tools, specifying the --type option - for example:
prefetch --type TenX --max-size 100000000 SRR5167880

Because the file was larger than the default limit of 20GB, I also had to increase it by specifying the --max-size option.
Alternatively, you can follow the official 10X Instructions to Download and Process BAM files of 1.3 Million Brain Cells.
I know this is an old thread but I'm hoping my answer will be helpful to anyone who has the same question (like me).

This is the simplest approach (if the process’s data is not available on GEO). Computing the UMI count matrix for 10x Genomics data is already implemented in Cell Ranger and if the same version is used it should reproduce the results of the study.

You should be able to parse out what you need using the tags in the .bam.  10xGenomics' website says what tags they add.
https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam
Going backwards would also involve parsing the tags, because you have to make two fastq files, and a simple bam -> fastq pipeline won't do that correctly.

Could you be more specific please? I'm very new to this procedure and bash scripting in general.

Yes, 10xGenomics single cell bams have the genes called in the bam tags.  https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam And no, samtools fastq will not make paired end fastqs suitable to rerun through cellranger, because the read data from the read containing the cell barcode and UMI will not have its own line in the bams.

@swbarnes2 There are indeed many tags, but isn't the UMI tag different than UMI counts? So how would I derive the expression matrix from these said tags? Thank you btw for your answer!

Each read has its UMI as a bam tag.  You can go through the bam line by line and figure out which reads were assigned to what cell barcode, UMI and gene.

This reverse process is also complicated by the fact that error-correction is performed on the barcodes and UMIs as well. For scRNA-seq, you typically start with 4 FASTQs: UMIs, barcodes, R1, and R2. There will likely be errors in those sequences, so barcodes/UMIs are corrected for by comparing to a whitelist of known input barcodes. You can take a look at the "Feature Barcoding Tags" section of the 10x support page @swbarnes2 listed to see what tags to look for. As mentioned, it's not a simple BAM -> FASTQ operation, but you should be able to do it

@James, that's not at all how the fastqs come out for 10xGenomics samples. (Or dropSeq either)  And feature barcoding is something new and extra; an old experiment probably won't have it.

My mistake, I was mixing up the scATAC-seq and scRNA-seq protocols. The scATAC protocol will have 3 or 4 FASTQs (https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/using/fastq-input), whereas scRNA has 2 or 3 (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/fastq-input). They both do barcode correction, but the input FASTQs are slightly different, and you're correct that the UMIs are probably only on newer datasets

@James UMIs have been used from at least cellranger 1.1.0 in 2016, if not earlier. Features are different, and were just introduced this year.

Cell Ranger bam files have the corrected and uncorrected barcodes and UMIs.  It is definitely possible to go from a Cell Ranger bam file to fastqs. 10X has a tool for doing it (https://support.10xgenomics.com/docs/bamtofastq)


text:
Identify non-coding regions from a genome annotation
I have this GTF file and I use the command below on a Linux machine to extract the coding regions of the genome:
awk '{if($3=="transcript" && $20=="\"protein_coding\";"){print $0}}' gencode.gtf

How I could do the inverse and keep only non coding regions?
I need non coding regions of human genome, but for asking my question here I referred to coding parts too.
After googling I see I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions.
Likely non-coding regions of genome are here:
https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.
UTR has its own feature in the gtf file. So you can do this:
$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

If the gtf file is compressed use this instead:
$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

BTW: Why are you using such an old release of gencode? The current version is v29.
As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with $ awk -v FS="\t" '$3=="CDS"' gencode.gtf. Reading the manual for bedtools complement is your task.

Sorry but your commands return nothing, I mean not working returning empty file

The gtf file the OP has linked to includes non-coding transcripts (LINCs, pseudogenes, tRNAs etc). I am guessing this is what they're after.

@FereshTeh if you run the command on the file you link to, it most certainly produces output. If it doesn't, you are using another file or have somehow modified the file.

This isn't a problem that's easily solved with awk. It's not like you're extracting a feature that's annotated in the GTF file. Instead, you want the empty space between annotated features.
A few years ago I wrote a program called LocusPocus for a similar task. It uses a gene annotation to break down a genome into gene loci and intergenic regions. It handles overlapping annotations and other weirdness pretty robustly. The output will include both coding regions and non-coding regions, but you can identify the intergenic spaces as those with iLocus_type equal to iiLocus or fiLocus.
Note: the --delta parameter will extend each gene/transcript by 500bp by default.
Caveat: the program only accepts GFF3 input by default. Hopefully it won't be too hard to convert your GTF to GFF3.
Another caveat: eventual interpretation of these data will depend on what features are annotated in the genome and which annotations you include vs ignore. Do you want your non-coding regions to include non-coding genes, or should these be treated separately? Some non-coding regions will be full of transposable elements and other repetitive DNA, while others will have enhancers, promoters, or other regulatory elements. It's important to tread carefully before you jump to any conclusions.

Absolutely brilliant name! :)

If you want all transcripts from that gtf file whose type isn't "protein_coding", you can use almost the same command, just change the == ("is") to != ("isn't"):
awk '{if($3=="transcript" && $20!="\"protein_coding\";"){print $0}}' gencode.gtf

Or, a simpler version:
awk '$3=="transcript" && $20!="\"protein_coding\";"' gencode.gtf

Note that this will not include any of the havana transcripts in the file, but I am assuming that's what you want since that's what your original command did.
Specifically, the command will return the following types of transcript (the numbers on the left are the number of such transcripts in the file):
awk '$3=="transcript" && $20!="\"protein_coding\";"{print $20}' gencode.gtf  | sort | uniq -c | sort -nk1
      1 "translated_processed_pseudogene";
      2 "Mt_rRNA";
      3 "IG_J_pseudogene";
      3 "TR_D_gene";
      4 "TR_J_pseudogene";
      5 "TR_C_gene";
     10 "IG_C_pseudogene";
     18 "IG_C_gene";
     18 "IG_J_gene";
     22 "Mt_tRNA";
     25 "3prime_overlapping_ncrna";
     27 "TR_V_pseudogene";
     37 "IG_D_gene";
     58 "non_stop_decay";
     59 "polymorphic_pseudogene";
     74 "TR_J_gene";
     97 "TR_V_gene";
    144 "IG_V_gene";
    182 "unitary_pseudogene";
    196 "IG_V_pseudogene";
    330 "sense_overlapping";
    387 "pseudogene";
    442 "transcribed_processed_pseudogene";
    531 "rRNA";
    802 "sense_intronic";
    860 "transcribed_unprocessed_pseudogene";
   1529 "snoRNA";
   1923 "snRNA";
   2050 "misc_RNA";
   2549 "unprocessed_pseudogene";
   3116 "miRNA";
   9710 "antisense";
  10623 "processed_pseudogene";
  11780 "lincRNA";
  13052 "nonsense_mediated_decay";
  25955 "retained_intron";
  28082 "processed_transcript";

You might also want to remove that "translated_processed_pseudogene" since that is actually translated into protein and is therefore technically coding:
awk '$3=="transcript" &&
     $20!="\"protein_coding\";" &&
     $20!="\"translated_processed_pseudogene\";"' gencode.gtf

Thanks a lot, really thank you for saving me I could not solve that myself.


text:
Separation of mixed plasmid DNA sequences post whole-plasmid sequencing
Imagine a DNA sample containing a mixture of different intact plasmids. These samples are sequenced using either MiSeq or HiSeq sequencing. Would it possible to assemble these plasmids post-sequencing as they would have been when sequenced individually?

Would you know the sequences? I mean, will you have reference sequences for each plasmid which you can use for assembling? How similar are the plasmid sequences? Will you have multiple ~100nt sequences that align perfectly to many different plasmids? Please [edit] your question and add more details.

Don't forget a template assembly, to confirm your results once you are clear what plasmids you are dealing with .. e.g. bowtie2, any paper on bacterial metagenomics will help you. Read length vs. sequence similarity is important

You would expect to have high coverage, given the plasmids are short, so de novo assembly would be likely very easy.  Given that each plasmid is present in different multiples, you would expect different coverage on each plasmid, so it might be best to approach it as a metagenome-type or transcriptome assembly, rather than a classic genome assembly.
Alternatively, Google search just told me there is a tool called plasmidSPADES designed for just your purpose.  SPADES itself is a reliable genome assembler, so I would imagine plasmidSPADES to be a good choice.  It does assume that there is also some whole bacterial genomes in the sequencing mix though.
As @terdon commented, it does depend on whether you have many identical copies of each plasmid, or many similar but not identical plasmids in the mix.

Yes I had the same ideas. I am very acquainted with SPADES,although never tried de novo assembly on a mixture of plasmids.

It looks as though plasmidSPADES expects to find some 'long' chromosomes in the mix, so it might be important to throw in some whole genome DNA into your mix of plasmids, to train the assembler.  You could do this after the experiment, by appending a fastq file of E.coli genomic DNA to your fastq file of plasmid sequences...

Coincidentally, the DNA sample is "contaminated" genomic DNA since it was extracted from an environmental sample. That should provide the necessary SPADES training?

If it's from an environmental sample might be better to try a metagenome assembler, maybe metaspades.  Alternatively, use kmer spectrum to try to pull the dataset apart based on abundance, and run the assembler on cleaned-up subsets of the reads

Let's take a step back and consider the "perfect" output for a de novo assembly algorithm. Ideally, you would like to see one complete sequence for molecule (chromosome, plasmid, etc.). In reality, this is difficult to achieve due to a couple factors.

By random chance some regions may have low coverage, meaning that there are an insufficient number of reads spanning the region to reconstruct it de novo.
A bigger problem is that some portions of the genome are repetitive and occur at many different locations, sometimes even on different chromosomes.

So in practice, it's very common that de novo assemblies can only recover partial fragments of each molecule, and sometimes fragments are erroneously fused due to shared repetitive content.
Back to your question: it's certainly possible to assemble distinct plasmids that were sequenced together without any special sample prep. How well this will work in practice will depend on the amount of coverage/redundancy in the sequenced reads and the amount of DNA shared between the plasmids. If they don't share much in common, and you sample to a sufficient depth of coverage, you shouldn't have any problem assembling them.

I think your comments are very informative so thank you so much!

@RoelofCoertze I'm glad it was helpful!


text:
Are duplicate variants against the VCF standard?
I have been trying to understand an error that the EBI's vcf_validator gives when run on my vcf file. Consider this minimal file:
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   0.00    .   .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

If I then run vcf_validator on it, it complains about:
$ cat $(LC_ALL="C" ~/vcf_validator_linux  -i test.vcf |& awk '/written to/{print $NF}' )
According to the VCF specification, the input file is not valid
Error: Duplicated variant chr12:529516:C>T found. This occurs 2 time(s), first time in line 6.

This is because the two variants are indeed the same. Changing
AACAC to AATAC at position 529514 just means change C to T at position 529516.
OK, but how is that against the VCF standard? What makes that an invalid VCF file? I am trying to automatically validate input VCF files and stop the analysis if they're invalid, but this seems like a perfectly valid file. Why is  vcf_validator not letting it through? The vcf-validator perl script from vcftools doesn't complain about this. So is it just a bug in  EBI's vcf_validator or is it really somehow illegal to have duplicate variants in a VCF file according to the specs?

What is the use case for having duplicate variant records in a VCF file?

@JohnMarshall none that I know of. But someone sent me a file with such duplicates and I was surprised that they caused the vcf_validator to consider the file invalid, so I asked if there was anything in the specs that would actually make such files invalid.

No. The only point the specification makes is that:

It is permitted to have multiple records with the same POS

I suspect this is intended to let the file list different variants occurring at the same position, but this incidentally allows duplicate variants.
This said, you probably don't want duplicate variants, because it will mess with any statistics you apply over the variants, and will unnecessarily clutter the file and slow down processing, if this happens a lot. You might also get variant annotations on one variant but not another, which can become a problem depending on how you query your variants.

Duplicate variants also point to an issue with upstream analysis. If the variants have different statistics (e.g. in this example one has qual 0.00, and the other has qual "."), then there is a bias in the counting. If the variants have the same statistics, then the lines are redundant.


text:
Mutation annotation to the heatmap
I have mutation data of patients.How do i add that mutation info or annotate those as column with my annotation information expression and mutation
To help myself i did try this Zuguang Gu but couldn;t able to do it, i would like similar to the one to show patient mutations .
any help or suggestion would be highly appreciated .
I wanted something like this , as here they annotated with batch ,type and cell clone. I would like to have something similar where i would like to see how mutation in patients affects expression profile .

df <- read.csv("HSC_pHSC_LSC_data/heatma_cluster.txt",header = TRUE,sep = ",")
        head(df)
        tail(df)
        dim(df)
        heat <- t(scale(t(df[,3:ncol(df)])))

        head(heat)

        rownames(heat) <- df$Gene
        ColAnn <- data.frame(colnames(heat))
        Sample = c(rep("HSC",4),rep("pHSC",12),rep("LSC",8))
        anno <- as.data.frame(Sample)
        categories <- data.frame(anno)
        colnames(anno) <- c("Condition")
        ColAnn <- HeatmapAnnotation(df=anno,which="col", col=list(Condition=c("HSC"="black","pHSC"="red","LSC"="dodgerblue4")))
        RowAnn <- data.frame(df$Family)
        colnames(RowAnn) <- c("Gene family")
        colours <- list("Gene family"=c("Cluster_one"="red","Cluster_two"="blue","Cluster_three"="brown","Cluster_four"="yellow",
                                       "Cluster_five"="green","Cluster_six"="maroon","Cluster_seven"="pink","Cluster_eight"="navyblue",
                                       "Cluster_nine"="darkgrey","Cluster_ten"="violet"))

        RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")
        boxAnnCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                            pch=".", size=unit(2, "mm"),
                                                            axis=TRUE,  axis_param = list(side = "left")),
                                       annotation_width=unit(c(1, 6.5), "cm"))

        boxAnnRow <- rowAnnotation(boxplot=row_anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),pch=".",
                                                            size=unit(3, "cm"), axis=TRUE, axis_param = list(side = "bottom")),
                                   annotation_width=unit(c(3), "cm"))
        pdf("ALL1_CLUSTER_DISEASE.pdf",width = 10, height = 20)
        #png("HSC_CMP_DIFF_HEATMAP.png", units="in", width=20, height=10, res=100)

        #dev.off()

        myCol <- colorRampPalette(c("navyblue", "white", "red"))(100)
        myBreaks <- seq(-2, 2, length.out=100)
        hmap <- Heatmap(heat,
                        name="Z-score",
                        col=colorRamp2(myBreaks, myCol),
                        heatmap_legend_param=list(color_bar="continuous",
                                                  legend_direction="vertical", legend_width=unit(5,"cm"),
                                                  title_position="topcenter", title_gp=gpar(fontsize=10, fontface="bold")),

                        #Split heatmap rows by gene family
                        split=df$Family,
                        #split = NULL,

                        #Row annotation configurations
                        cluster_rows=TRUE,
                        show_row_dend=FALSE,
                        row_dend_side = "left",
                        #row_title="Transcript", #overridden by 'split' it seems
                        row_title_side="left",
                        row_title_gp=gpar(fontsize=8,fontface="bold"),
                        row_names_gp = gpar(fontsize = 8),
                        show_row_names=TRUE,
                        row_names_side="right",
                        row_title_rot=0,

                        #Column annotation configuratiions
                        cluster_columns=TRUE,
                        show_column_dend=TRUE,
                        column_title="Samples",
                        column_title_side="top",
                        column_title_gp=gpar(fontsize=10, fontface="bold"),
                        column_title_rot=0,
                        column_names_gp = gpar(fontsize = 20, fontface="bold"),
                        show_column_names=TRUE,

                        #Dendrogram configurations: columns
                        clustering_distance_columns="euclidean",
                        clustering_method_columns="complete",
                        column_dend_height=unit(10,"mm"),

                        #Dendrogram configurations: rows
                        clustering_distance_rows="euclidean",
                        clustering_method_rows="complete",
                        row_dend_width=unit(15,"mm"),

                        #top_annotation_height=unit(0.5,"cm"),
                        top_annotation=ColAnn,

                        #bottom_annotation_height=unit(2, "cm"),
                        bottom_annotation=boxAnnCol)
        hmap + RowAnn

        dev.off()

The image i get after running the code is something like this sample numbers vary, as i have only HSC ,pHSC and LSC samples

Hi, could you please [edit] the question to fix the punctuation (commas before space, n't instead of n;t and move i to I) and the code indentation? It will be easier to read. Also could you please add the final image and a description of what do you want to change? Without seeing it and having your data is hard to guess how we can help

okay will do it

What exactly are you trying to do? Still don't understand..

i have expression data and mutation data , i wanted to add the mutaation annotation to for each patient in heatmap as well ..

edited the question..i hope now it's bit more clear

Could you please add the image you get after all this code?

okay...i will .its a normal heatmap using complex heatmap what i get..i will be posting it..

I assume that the x-axis is patient and the Y-axis is the gene?  The results sort of look clear, but threading together patients and "batch type cell clone" as well as "HSC ,pHSC and LSC" is a mystery to me.  I assume someone here understands this.

no i don;t have batch type cell close, instead i have mutation data ..yes x and y axis are patient and gene . HSC is stem cell ,pHSC is patient stem cell and LSC is leukemic stem cell.

You are trying to "map" mutation data to your heatmap, however, there is no mutation data in your code as far as I can see! You specify sample names within a data frame called anno but this data frame has only one column, Condition. All you have to do is add another column, mutation, to the anno data frame. Moreover, you have to pay attention to the order of your data and your meta data. The rownames of the meta data (anno data frame in your case) should be the same as the colnames of your gene expression data frame.
I also see that you create a row annotation object RowAnn, which is not used in the Heatmap() call. Similarly, there is a categories data frame which is not used either. Your code would be easier to read without these unused variables.

i did figure it out after spending some time on complex heat-map documentation


text:
awk for motif splicing from an amino acid alignment
I'm looking to splice out specific amino acid alignment positions as a continuous motif from an alignment (~15 000 sequences). The code is written in Perl (I agree not many developers now) but splicing method is awk which a lot of you know. Please note the hash is "tied" (not shown), this just means the order of the sequences is preserved.  I can get the awk code below to return single amino acid positions.
I am struggling to get a loop working in awk.
I agree the Perl wrapper makes it more complicated, but its part of a huge data pipeline (all written in Perl).
#!/usr/local/ActivePerl-5.20/bin/env perl
#use strict;
use warnings;

# I loaded the alignment into a hash, where key = sequence and value = protein
my %dna = (
    virus1 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF",
    virus2 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF",
    virus3 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF"
);

# I define the start and end of the alignment positions I am splicing
my $n = 1;
my $z = 5;

# I turn the alignment into a string with separated by carriage returns
my @array = map {"$dna{$_}"} keys %dna;
my $string = join "\n", @array;

# This is the awk loop, note in Perl the "``" means issue a shell command

my $splice = `echo "$string" | awk 'BEGIN {FS=""} {x = $n; while (x < $z); print \$ x; ++x}' 2>1`;
print $splice;
__DATA__

Expected output
LGVLV
LGVLV
LGVLV

Observed output
I get the original sequences returned
I am familiar with Python/Perl solutions of line by line processing but awk is more efficient IMO.

Could you [edit] and show the output you're expecting? I speak both Perl and awk, but I don't get what output you want. Also, I really doubt that calling awk will be more efficient than doing it natively in perl. You don't seem to be using anything that awk is particularly good at.

Amended, thanks. I will post the solution I got working, but this is because I don't understand awk

The main problem is that your while only encloses the first statement after it, and not the second. You need { } for that. If you could explain what that awk is supposed to be doing (I'm pretty sure you've copied it wrong, is there really a space between \$ and x?), I should be able to post an answer. Do you just want the first 5 characters of each sequence? Don't you want positions $n to $z (3-7)?

I think you are trying to use awk to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:
my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

That, when added to your perl script, produces:
$ foo.pl
VLVIL
VLVIL
VLVIL

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.
Now, the reason your version didn't work, is that the while was like so:
while (x < $z); print \$ x; ++x}

Assuming the \$ x was a typo, the problem is that the ++x isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need {}:
while (x < $z){print \$x; ++x}

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use awk's substr command:
  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

So, in your script, that would be:
my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use awk instead. It will not be faster or more efficient than the same thing in perl:
print substr($_,$n,5) . "\n" for @array;

By calling awk you are also calling the shell which will be slower by definition, especially if you're echoing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the awk loop I gave above, or the awk+substr, or the native perl substr:

awk loop, average time taken to run the same thing 100 times: 0.0525
awk substr, average time taken to run the same thing 100 times: 0.0522s
perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:
print substr($dna{$_},$n,5) . "\n" for keys(%dna);

That averaged at 0.0376s over 100 runs.

Huge thanks @terdon. You are right I should have performed a runtime check :-( such as use Timer::Runtime. I worked on the implicit assumption awk was faster for this operation. David Cross "lectured" on runtime in detail in "Data Munging for Perl": There is no denying that 0.0525 vs. 0.0376s is 30% faster.

@MichaelG. in my experience, awk is blindingly fast when your data are organized in fields because of its implicit splitting. If not, however, it doesn't give any real advantage over a language like Perl which is simply brilliant for text parsing.

If you have the sequences in a seperate file like this:
$ cat input
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF

you can use a simple cut to fetch a substring:
$ cut -c1-5 input
LGVLV
LGVLV
LGVLV


text:
CollectHsMetrics base coverage output has overlapping targets
I am trying to get base coverage information, and am using the --PER_BASE_COVERAGE output from Picard tools to get a text file of base level coverage. However, some of my targets are being grouped together. In the targets file, these are separated by one base as such:
12      12258595        12258695        +       TARGET1
12      12258696        12258795        +       TARGET2
12      12258796        12259791        +       TARGET3

However, when I look at my covg file, I see the following:
chr    pos    target   covg
12  12258595        TARGET1|TARGET2|TARGET3 510

Why would HSmetrics group these targets together, and how can I separate them?

Question is also asked on biostars.

Hello @Blaze9, I don't understand your question. Would you maybe explain bit what are the targets? Are they read groups? (disclaimer, I am not a Picard tools user); I usually use samtools for this.

Hi @KamilSJaron, the targets are either tiling regions or exons/introns. I can checkout samtools to see if it provides better results. The per_target_coverage output does not combine the targets which is the weird part. Thanks!

Ah, I see. Maybe you want to [edit] your question and add as many details as you can (such those you wrote in the comment) :-).


text:
Why can I not install snakemake on my SLURM computer even though I can find it in the bioconda channel?
I searched for the package 'snakemake' on my SLURM cluster using:
conda search --channel bioconda snakemake
and I get many versions, up to 5.4.2.
I then try to install it using:
conda install --channel bioconda --yes snakemake=5.4.2
but it fails with the following message:
PackagesNotFoundError: The following packages are not available from current channels:

  - snakemake

Current channels:

  - https://repo.anaconda.com/pkgs/main/linux-64
  - https://repo.anaconda.com/pkgs/main/noarch
  - https://repo.anaconda.com/pkgs/free/linux-64
  - https://repo.anaconda.com/pkgs/free/noarch
  - https://repo.anaconda.com/pkgs/r/linux-64
  - https://repo.anaconda.com/pkgs/r/noarch
  - https://repo.anaconda.com/pkgs/pro/linux-64
  - https://repo.anaconda.com/pkgs/pro/noarch

Why?
Why is the bioconda channel not included in the current channels? I thought that was what the --channel bioconda did in my install command.
Note that conda install --channel bioconda --yes samtools=1.9 worked fine and installed samtools successfully.

I've heard some folks warn that conda is sensitive to the order in which channels are set/configured. I'm not exactly sure what the details are, but maybe that gives you somewhere to start troubleshooting?

I created a file .condarc in my home directory with the following content:
channels:
  - conda-forge
  - bioconda
  - defaults

and it works now. I don't understand why this is necessary though. I thought the channels in the .condarc file would be used if I did not specify them in my command.

This seems to be a bug in the version of conda you're using, it works in 4.6.4. Note, however, that you needed to add -c conda-forge to your original command anyway, since many of the snakemake dependencies come from there.

good point about the dependencies. I didn't know about that. I just always give the channel of the package I want to install but it sounds like I also need to give the channel of dependencies. So I'll keep using -c conda-forge -c bioconda -c r in all my conda install commands.


text:
Are there any Genome in a Bottle-like resources for non-humans (especially for invertebrates)?
Genome in a Bottle is an excellent resource that provides many types of DNA and RNA sequencing reads for a single individual/cell line to test genome assembly and analysis tools. For example there are Illumina WGS reads, Oxford Nanopore, PacBio, and 10X reads available for download at high coverage depth. In addition, there are pre-aligned bam files and verified variant calls.
Is there any aggregation of similar robust single-individual sequencing resources for non-human, non-mammalian metazoan species? The only way I can consider finding such resources at the moment is going through SRA projects one at a time and hoping for the best.

Have you checked sites like http://www.flybase.org/ and https://wormbase.org/?

In addition I would like to add "Vectorbase", the curator "Dan" is a good guy (I forget his surname right now) and it's a good team. However, please read my rather somber assessment of arthropod genomics first.

We appear to be talking about lower metazoans, e.g. not protists, or yeasts (the latter I suspect are really good for your question). I don't know anything about worms.
Its important to understand that arthropods are several (soccer) divisions lower than mammals in terms of the question you are asking, except for two model organisms.
Flybase, Drosophila melanogaster is the most important species for your purposes. The second most important species within the Class Insecta will be Anopheles gambiae (malaria mosquito).
Personally I would stick to the main players because overall whilst there are a lot of insect genomes (because there are a lot of insects), their genomes remain rough because of the difficulties in assembly and requirement to physical mapping their genomes in comparison to their genome size. Also except for D. melenogaster the (incredible) tardiness in producing a draft genome. For example the idea of PacBio is in its infancy for many entomologists: and this would help them a lot in assembling their genomes correctly. I agree things have picked a lot recently, but I remain a skeptic.
Within the arthropods I would be very surprised if the crustacea, myriapods or Chelicerata (e.g. spiders) were in better condition than hexapods.
If you are interested in protists certainly Plasmodium falciparum (African malaria) has the sort of information you seek (Plasmodb).

Summary
I don't want to discourage you, insects are great, but you need to be selective if you are from a human/mammalian genetics background.

Does FlyBase or Plasmodb actually have the sort of data GIAB offers? High quality reads files and validated variants for those reads? And there are draft (and not so draft) genomes for dozens of insects by now. I was involved in the 12 Drosophilas paper, more than a decade ago, and have worked with A. gambiae, A. aegypti, T. castaneum, A. melifera among others. I've since left the field of comparative genomics, but I am surprised at your claim of a scarcity of genomes given there were already several available a decade ago. Have they since stopped working on them?

Wait, wait ... . I'm not claiming scarcity of genome - no way - (there are lots of insect genomes), I am claiming quality of genome will still be absent. Drosophila is an exception, with D. melanogaster exceptional quality. For example, physical mapping against the polytene chromosomes was certainly lacking for most insects but is easily doable and would really help assembly. I was in on tsetse and a hemipteran genome.  GIAB I don't keep up but I'd be surprised if melanogaster was not up to this standard.

Ah, OK, that makes more sense. It was the "tardiness in producing a draft genome" that threw me. As for GIAB, it isn't so much the quality as the type of data. They offer raw reads and.or aligned bam files and, most importantly, a "truth set" of variants that should be found in those reads. So it's a great (well, it's the only one, anyway) tool for validating variant calling pipelines. I don't know if such a dataset is available for non-human organisms.

@MichaelG. I work with marine invertebrate metazoans! Thanks for your answer.


text:
Aligning nucleotide sequences in APE software
I am very new to APE software. I am trying to align complementary oligo-DNA strands using APE on macOS. I have both the forward and reverse sequence of an oligo-DNA insert in two separate files. But I am unable to align them as complementary sequences and show the sticky ends.
Is it possible with APE? I tried Align two sequences option. But couldn't get the expected result.

I added a link to the tool I am guessing you are using. Please correct it if it's the wrong tool. Also, please show us the result you get so we understand the problem. Does the program crash? Does it not align? Does it give a bad alignment? Does it just not give the one you were expecting?

The link is important because of confusion with R's "ape" package. I'll check it out. Are you using Mac or PC?

@Michael G. I am using Mac

Thanks @Vikram . Yes, me too, I'll take a look.

Could you upload the files you are trying to use as input? Or paste the text in?

@Science123, Don't know if this is what you are looking for but here is the code for a minimal working example of using APE to produce a tree from a multiplesequencealignment loaded from a file; https://bioinformatics.stackexchange.com/a/14160/9311


text:
how to obtain GRCh38 bed and symbol files for HOCOMOCOv11
I am trying to a recently published software called diffTF and contacted the authors on how to obtain the input HOCOMOCO files for GRCh38. While I wait for the answer, I thought I would post the question here as well:
The tool requires bed files for each of the transcription factors (TF) in the genome, such as this bed file (TAB separated), which seems to give the locations of different TF binding sites in the genome:
chr1    3457132     3457141     GTTGCGTGA   1266   +
chr1    4595374     4595383     GTTGCGTGA   1266   -
chr1    4996631     4996640     TTTGCGTGC   1276   +
chr1    5005357     5005366     TTTGCGTGC   1276   +
chr1    5170456     5170465     TTTGCGTGC   1276   -
chr1    5567302     5567311     GTTGCGTGA   1266   -
chr1    8676279     8676288     TTTGCGTGC   1276   -
chr1    8859719     8859728     GTTGCGTGA   1266   +
chr1    9600920     9600929     GTTGCGTGA   1266   +
chr1    9650971     9650980     GTTGCGTGA   1266   +
...

Then there is a symbols file to map the TF to their gene(?) that looks like this:
SYMBOL  ENSEMBL HOCOID
Ahr     ENSMUSG00000019256      AHR
Aire    ENSMUSG00000000731      AIRE
Alx1    ENSMUSG00000036602      ALX1
Ar      ENSMUSG00000046532      ANDR
Tfap2a  ENSMUSG00000021359      AP2A
Tfap2b  ENSMUSG00000025927      AP2B
Tfap2c  ENSMUSG00000028640      AP2C
Tfap2d  ENSMUSG00000042596      AP2D

I googled for HOCOMOCO and found a bunch of pwn or words files, but none of them seems to be similar to the bed file above. Any ideas where to get these for GRCh38?

What have you tried besides contacting the authors? I am not familiar with this type of files, but it doesn't seem to be too strange, so perhaps they assumed it is fairly easy to download and transform from other files into these...

I guess you're right that you need predictions of transcription factor binding site (TFBS) through genome.
There're two options to get them:

You can generate predictions using any motif scanning tool. E.g. you can use our tool SPRY-SARUS which is able to generate bed-file in necessary format.
In this case, you can choose more or less strict thresholds for each prediction (i.e. more false positives vs more false negatives). Standard thresholds for HOCOMOCO motifs can be found at its downloads page (we, as HOCOMOCO authors, recommend thresholds corresponding to P-value of 0.0005).
I had no experience with diffTF but it seems that authors prepared TFBS predictions and put them on their site. Please, take a look at this page of manual. If you don't need extra-flexibility in the generation of TFBS predictions, it can be the simplest way to proceed.