yctuga
commented
3 years ago BINF8211 HW3
Yun-Ching Tsai
Q.0. shell script cufflinks -p 8 -g $gtf -b $genome -u -o $result $result/PE_2_hg38_sorted.bam
-p: number of threads used for analysis -g: the reference annotation file for assembly -b: run with new bias detection and correction algorithm which can significantly improve accuracy of transcript abundance estimates -u: do an initial estimation procedure to more accurately weight reads mapping to multiple locations in the genome -o: output file destination
http://cole-trapnell-lab.github.io/cufflinks/cufflinks/index.html
Q.1. How many genes, transcripts and exons are identified by Cufflinks?
genes: $ cat genes.fpkm_tracking |cut -f1 |wc -l 57083
$ cat genes.fpkm_tracking |cut -f1 |grep -c ^CUFF 18185 (genes assembled by CUFFLINKS)
transcripts: $ cat isoforms.fpkm_tracking |cut -f4 |wc -c 2805187
$ cat isoforms.fpkm_tracking |cut -f4 |grep -c ^EN 60444 (omit CUFFLINKS)
exon: $ cat transcripts.gtf |cut -f-1,3,4,5,7 |awk 'BEGIN{FS="\t";OFS="\t"}{if($2=="exon"){print $0}}' |sed 's/\t/**/g' |wc -l 1493389
Q.2. What is the highest expressed gene(s) and transcript(s)?
genes: $ cat genes.fpkm_tracking |cut -f1,10 |sort -rn -k2,2 |more
CUFF.17999 4237.89 (maybe not correct) ENSG00000160182 877.525
transcript(s): $ cat isoforms.fpkm_tracking |cut -f1,10 |sort -rn -k2,2 |more
CUFF.17999.3 4237.89 ENST00000293308 4006.55
!/bin/bash
SBATCH --job-name=hw3_Cufflinks
SBATCH --partition=batch
SBATCH --ntasks=1
SBATCH --cpus-per-task=8
SBATCH --time=24:00:00
SBATCH --mem=20gb
SBATCH --output=hw3.%j
SBATCH --mail-user=yct@uga.edu
SBATCH --mail-type=ALL
define directory
data='/home/yct/BINF8211/HW2/data'
sam_input='/home/yct/BINF8211/HW2/result' result='/home/yct/BINF8211/HW3/result' genome='/work/binf8211/instructor_data/HomeWork2/source/hg38.fa' gtf='/work/binf8211/instructor_data/Homo_sapiens.GRCh38.102.chrV1.gtf'
module load
module load HISAT2/2.1.0-foss-2019b
module load SAMtools/1.10-GCC-8.3.0 module load Cufflinks/2.2.1-foss-2019b
sort bam file
samtools view -bS $sam_input/PE_2_hg38.sam > $result/PE_2_hg38.bam
samtools sort $result/PE_2_hg38.bam $result/PE_2_hg38_sorted
mkdir $result/tmp samtools sort -T $result/tmp/ -o $result/PE_2_hg38_sorted.bam $result/PE_2_hg38.bam rm -r $result/tmp
cufflinks
cufflinks -p 8 -g $gtf -b $genome -u -o $result $result/PE_2_hg38_sorted.bam