How to prepare .gct file for GSEA software?

Question

[zekunmu]

Hello~

I am working on my NGS data for gene set enrichment analysis now. I wonder how can I convert my .fastq CleanData received from sequencing company into .gct file for GSEA software?

I know some software like GenePattern can convert some file format such as .arff and .pcl into .gct file, but it seems no such tool is available for .fastq file!

I'd appreciate it a lot if some could help!!!

[moxu]

As WouterDeCoster said, you are really skipping too many steps.

If all you care about is GSEA, you can rank your genes using some criteria (such as from DEG analysis) and prepare a rank file. You can then feed the ranked gene list to GSEA.

Answer 1

[Barry Digby]

For future users who end up here, run this script:

bash format_GCT.sh norm_counts.txt output_prefix

Script assumed you have saved normalised counts from DESeq2 in R:

norm_counts <- counts(dds, normalized=TRUE)
write.table(norm_counts, file="norm_counts.txt", append = TRUE, sep="\t", row.names = TRUE, col.names = TRUE, quote = FALSE)

Here is the format_GCT.sh bash script:

#!/usr/bin/env bash

## Script to convert DESeq2 normalized counts to GCT format

if [ "$1" == "-h" ]; then
  echo "Usage: bash `basename $0` [FILE] [OUT_PREFIX]"
  echo
  echo "Example: bash `basename $0` norm_counts.txt GSE118959"
  exit 0
fi

if [ "$1" != "" ]; then
    echo "Counts File: $1"
else
    echo "Counts File is missing"
fi

if [ "$2" != "" ]; then
    echo "Output Prefix: $2"
else
    echo "Output Prefix is missing"
fi

FILE=$1
OUT=$2

TMP=`wc -l $FILE | awk '{print $1}'`
LEN=`expr $TMP - 1`
TMP=`awk '{print NF}' $FILE | sort -nu | tail -n 1`
SAMPLES=`expr $TMP - 1`


cat $FILE \
| awk -F"\t" '{if(NR==1) $1="NAME"FS$1}1' OFS="\t" \
| awk '{$1 = $1 OFS (NR==1?"Description":"na")}1' \
| sed 's/ /\t/g' \
| awk -v LEN=$LEN -v SAMPLES=$SAMPLES \
 'BEGIN{print "#1.2""\n"LEN"\t"SAMPLES}1' > ${OUT}.gct

[tianqi.ma1994]

Thanks! It's really helpful. I tried this script and it works but the row names and column names in the output file contains quote for each gene/cell. How to remove the quotes as GSEA gene set doesn't use quote.

[Barry Digby]

When saving the normalized counts in R set quote = FALSE .

Or just remove them using sed sed 's/"//g' your_file.gct > no_quotes.gct

Answer 2

[WouterDeCoster]

You are skipping an awful lot of steps of your analysis. Your fastq data first needs to be aligned to the genome, reads need to be counted and you need to perform differential expression analysis. The results of the last step can be used in GSEA.

[zekunmu]

Thank you so much! I have worked it out, what you said is really true. It was my first time working on NGS data. I finally found out the sequencing company actually processed the data somewhat, such as alignment and reads counting.