htseq-counts output merge into one matrix ??

I just need a little help to merge my all features counts into one matrix. I have counted features using htseq-counts and now want to merge into one file like....

ID c1 c2 c3..............t1 t2 t2.......etc

The problem is that i have 36 sample and corresponding counts so i am lazy to write ...

paste *.txt | cut -d -f1,2,4,6,8...............................72 (because each file contain two column i.e ID and counts )

and moreover, apart from laziness, one problem also happened here that while pasting its not in sequence like :

paste 1.txt 2.txt 3.txt............

Help me please...

rna-seq gene expression counts bioinformatics • 15k views

ADD COMMENT • link •

modified 8 months ago by Alex Nesmelov • 190 • written 6.4 years ago by bioinforupesh2009.au • 130

Finally got a solution to this problem

#!/bin/bash

FILES=$(ls -t -v *.txt | tr '\n' ' ');

awk 'NF > 1{ a[$1] = a[$1]"\t"$2} END {for( i in a ) print i a[i]}' $FILES > merged.tmp

This woks like you asked.

ADD REPLY • link written 6.4 years ago by EagleEye ♦ 6.8k

Thank you man for your kind reply but i have already done, i am sorry for the delay. Was just needed to change the file name like 1.txt -> 01.txt, 2.txt -> 02.txt....................... etc and after i used your suggested awk code and worked fine.

By the way thanks a million.

ADD REPLY • link modified 6.4 years ago • written 6.4 years ago by bioinforupesh2009.au • 130

6.4 years ago by

EagleEye ♦ 6.8k

Sweden

EagleEye ♦ 6.8k wrote:

This could be the simple way for as many as files:

awk 'NF > 1{ a[$1] = a[$1]"\t"$2} END {for( i in a ) print i a[i]}' read_counts/*.gene_counts.htseq > merged.htseq

ADD COMMENT • link written 6.4 years ago by EagleEye ♦ 6.8k

This outputs: Merge htseq files into one based on geneNames. Sample output

ENSG00000251243.1   0   0   0   0   0
ENSG00000240637.2   0   0   0   0   0
ENSG00000227366.1   0   0   0   0   0
ENSG00000104903.4   665   216   884   3254   689
ENSG00000228626.1   0   0   0   0   1

Forgot to mention, the order of pasted column will be in FileName order.

Like if you do listing using unix command 'ls -l read_counts/'

You will get the order in which it has merged.

ADD REPLY • link modified 6.4 years ago • written 6.4 years ago by EagleEye ♦ 6.8k

Dear Santhilal,

Your code is working as mine too...but the main problem what i mentioned on question is.......ITS NOT IN SEQUENCE.

I NEED...

ID C1 C1 C3.............UPTO 36

your code or mine work as terminal shows the directory by hit of ls

10.txt

11.txt

31.txt

12.txt

13.txt

4.txt

3.txt

1.txt

20.txt

...

....

36.txt

However, i tried by

for i in $(seq 1 36); do awk 'NF > 1{ a[$1] = a[$1]"\t"$2} END {for( i in a ) print i a[i]}' $i.txt > merged.htseq; done

but its not working ???

ADD REPLY • link modified 6.4 years ago • written 6.4 years ago by bioinforupesh2009.au • 130

you can try this

ls -l -v | awk '{print $9}'

and then AWK. I don't know whether it will work. But you can try it.

ADD REPLY • link written 6.4 years ago by EagleEye ♦ 6.8k

Thanks but while using this only its work but not with advised final awk code.

i used...but not worked for me.

ls -l -v | awk '{print $9}' | awk 'NF > 1{ a[$1] = a[$1]"\t"$2} END {for( i in a ) print i a[i]}' *.txt > merged.htseq

however, i have done by my way........

Thanks anyway....very kind of you. :))

by the way.....

Wishing you a very happy and prosperous Diwali!!

ADD REPLY • link written 6.4 years ago by bioinforupesh2009.au • 130

Could please provide your working solution?

ADD REPLY • link written 3.5 years ago by Ric • 330

6.4 years ago by

Devon Ryan ♦ 98k

Freiburg, Germany

Devon Ryan ♦ 98k wrote:

I wrote this a while ago. It can also change the file labels.

Note that it's so easy to just do this in R that there's really little reason to manually merge files like this.

ADD COMMENT • link modified 6.4 years ago • written 6.4 years ago by Devon Ryan ♦ 98k

and if you use DESeq2 you can simply give the file names in input ( see DESeq2 vignette 1.2.4 HTseq input )

ADD REPLY • link written 6.4 years ago by Nicolas Rosewick ♦ 9.3k

Exactly, and even if you're using something like edgeR, it's simple enough to just use the same method employed by DESeq2 (namely, lapply() reading files and then lapply a function to get the counts out).

ADD REPLY • link written 6.4 years ago by Devon Ryan ♦ 98k

i prefer first by terminal or bash.

Thanks by the for your valuable answer.

ADD REPLY • link written 6.4 years ago by bioinforupesh2009.au • 130

2.7 years ago by

Arindam Ghosh • 340

Finland

Arindam Ghosh • 340 wrote:

I managed to do this by R. The script is available here. It works for me. Please try and check if it works fine.

ADD COMMENT • link written 2.7 years ago by Arindam Ghosh • 340

4.6 years ago by

bmahesh07 • 0

United States

bmahesh07 • 0 wrote:

I found a useful R script (htseq-combine_all.R script) for merging HT-seq counts at the following resource http://wiki.bits.vib.be/index.php/NGS_RNASeq_DE_Exercise.4

ADD COMMENT • link written 4.6 years ago by bmahesh07 • 0

I tried the the above R script (htseq-combine_all.R script) but I ran into an error

ADD REPLY • link written 3.5 years ago by Ric • 330

You can use my R script on https://github.com/AAlhendi1707/htseq-merge/blob/master/htseq-merge_all.R

ADD REPLY • link modified 23 months ago • written 23 months ago by Ahmed Alhendi • 50

2.1 years ago by

oriolebaltimore • 160

United States

oriolebaltimore • 160 wrote:

Is it considered best practice to sum up reads across replicates or to run HTSeq on merged BAM File after ccombing fastqs. Thanks Adrian

ADD COMMENT • link written 2.1 years ago by oriolebaltimore • 160

No, you should use hts-seq on each replicates otherwise there is no meaning of replicates. By the way, you could use "featureCounts tool" for each bam at a time. Then, further you don't need to merge as its give you a matrix of count for each sample of bam.

Hope this helps Good luck.

ADD REPLY • link written 2.1 years ago by unique379 • 90

23 months ago by

Ahmed Alhendi • 50

University of Southampton, UK

Ahmed Alhendi • 50 wrote:

This is easily done using R script! You can use my R script on https://github.com/AAlhendi1707/htseq-merge/blob/master/htseq-merge_all.R

For more details about how to generate and combine HTSeq outputs into one read count matrix in R https://ahmedalhendi0.wordpress.com/2019/03/20/combine-htseq-outputs-into-one-read-count-matrix-in-r/

ADD COMMENT • link modified 23 months ago • written 23 months ago by Ahmed Alhendi • 50

Hello, I followed your script and ran it on Linux, but I got an error message. My input files were in the same directory where I run R script, but it said cant read files. How can I fix it?

Thx .

ADD REPLY • link written 6 months ago by Aynur • 40

Recheck on path and file format.

ADD REPLY • link written 6 months ago by Arindam Ghosh • 340

Hello Dr.Alhendi, I am totally new to this bioinformatics, so please bear with my stupidity here, if any. Thanks. I did check again. Here is how I used the Rscript. 1. I could not use it on Mac OS terminal. 2. I cant use it on Mac OS R. 3. I tried to use it on Linux cluster as follows.

`Rscript htseq-merge_all.R /kuacc/users/akashgari19/HTSeq  HTSeq_Merged.txt`

R script is located on the same directory. My HTSeq out files are located on this directory:

`/kuacc/users/akashgari19/HTSeq`

My files are named below : CON-1_HTSeq.txt treatment-1_HTSeq.txt ....... total 10 files.

`head CON-1_HTSeq.txt.                                                                     ENSMUSG00000000001   3141 
 ENSMUSG00000000003 0.                                                                  ENSMUSG00000000028  854

ENSMUSG00000000031 822918 ENSMUSG00000000037 1 ENSMUSG00000000049 3 ENSMUSG00000000056 3897 ENSMUSG00000000058 1399 ENSMUSG00000000078 1766 ENSMUSG00000000085 706` When I run the script I got the following errors. [1] "############################" [1] "## Files to be merged are: ##" "CON-1_HTSeq.txt" "CON-2_HTSeq.txt" ........

file read START ######### Aug_12_2020_13_18_31_+03" Error in file(file, "rt") : cannot open the connection Calls: read.table -> file In addition: Warning message: In file(file, "rt") : cannot open file '/kuacc/users/akashgari19/HTSeq/CON-1_HTSeq.cnt': No such file or directory. Execution halted

Yes, I do not have files with .cnt extension.
Please help me, how can I solve this. Thanks.

ADD REPLY • link modified 6 months ago • written 6 months ago by Aynur • 40

8 months ago by

Alex Nesmelov • 190

Alex Nesmelov • 190 wrote:

I decided to revive the thread adding my short tidiverse-based solution in R. Assuming we are in R session in the directory with htseq-counts files:

library(tidyverse)
counts_files_pattern = "\\.counts\\.txt$"
counts_table <-
  list.files(pattern = counts_files_pattern) %>%
  #apply an import function to the listed files
  map(function(x) { 
        print(x)
        read_table2(x,
                #rename second column in accordance with the file name
                col_names = c("ID", x))
      }
      ) %>%
  #merge all tables into one
  reduce(left_join, by = "ID") %>% 
  # clean up column names
  rename_all(gsub,
             pattern = counts_files_pattern,
             replacement = '') %>%
  #move ID column to row names
  column_to_rownames("ID")

ADD COMMENT • link written 8 months ago by Alex Nesmelov • 190

Hello, I used this code on R , but got an error. I ran the script on the same directory I had my count files. Got the following error. Error: .x is empty, and no .init supplied.

Can you help me? Thanks.

ADD REPLY • link written 6 months ago by Aynur • 40

Please log in to add an answer.

Similar posts • Search »