二代测序数据处理（二）从CleanData到表达矩阵

第一步 准备好数据

• CleanData：由测序公司提供

内有 name_R1.fq.gz 和 name_R2.fq.gz

• 对应物种的索引文件：二代测序数据处理（一）数据格式说明

大鼠的索引文件

• 对应物种的染色体模板

上面的染色体模板命名不规范，使用下面脚本修改

#!/bin/bash
for file in `ls /home/lxb/gene_template/rat/gene/Rattus_norvegicus.Rnor_6.0.dna.chromosome.*.fa`
do
    newFile='/home/gene/xuchen/ratgene/'`echo $file  sed 's/\(.*\)\/.*\.\([0-9XYMT]\+.fa\)$/\2/'`
    echo $newFile
    cp $file $newFile
done

规范后的染色体命名

• 对应物种的GFT注释文件

两份注释文件的差别，一般使用不带chr的文件

第二步 生成排好序的bam文件

• nohup /home/gene/xuchen/gp1.bash > /home/gene/xuchen/1.log 2>&1 &

#!/bin/sh
#设置CleanData存放目录
export CLEAN=/home/gene/upload/xc-lz/gp
#设置idx存放目录
export HISAT_IDX=/home/lxb/gene_template/rat/idx/rat.hisat2.idx
#设置这一步的输出目录
export WORK=/home/gene/xuchen/output_gp
  
export CDIR=$(basename `pwd`)
echo $CDIR
echo $CLEAN
for file in $CLEAN/*
do
echo $file
SAMPLE=${file##*/}
echo $SAMPLE
r1=${SAMPLE}"_R1.fq.gz"
r2=${SAMPLE}"_R2.fq.gz"
echo $r1
echo $r2
hisat2 -p 12 --dta-cufflinks -x $HISAT_IDX -1 $CLEAN/$SAMPLE/$r1 -2 $CLEAN/$SAMPLE/$r2 -S $WORK/$SAMPLE.sam
samtools view -bS $WORK/${SAMPLE}".sam" > $WORK/${SAMPLE}".bam"
# 内存使用量为 6000M * 6 = 36G
samtools sort -m 6000M -@ 6 -o $WORK/${SAMPLE}".sorted.bam" $WORK/${SAMPLE}".bam"
#没用的东西删掉
rm $WORK/${SAMPLE}".sam"
rm $WORK/${SAMPLE}".bam"
done

• 查看 1.log 确保匹配率高于80%，下图为 95.01%

• 最终结果

第三步 生成基因矩阵

• nohup /home/gene/xuchen/gp2.bash > /home/gene/xuchen/2.log 2>&1 &

#上一步的输出目录
export CLEAN_OUT=/home/gene/xuchen/output_gp
cd $CLEAN_OUT
#规范命名后的染色体模板数据
export GRCh38=/home/gene/xuchen/ratgene
#GTF数据，一般使用不带chr的
export GTF=/home/lxb/gene_template/rat/gtf/Rattus_norvegicus.Rnor_6.0.90.gtf
#整理样本列表, 空格分隔
export SAMPLES="A34 C15 C41 C8 M7 X24 X28"
#用来生csv的表头，逗号分隔，与样本列表对应
export LABLES="A34,C15,C41,C8,M7,X24,X28"
#设置输出目录
export WORK=/home/gene/xuchen/diffout_gp

export LIST=""
for file in $SAMPLES
do
LIST=$LIST" "$CLEAN_OUT/$file".sorted.bam"
done
echo $LIST

cuffdiff -L $LABLES -o $WORK -b $GRCh38 -p 12 -u $GTF $LIST

• 查看处理进度 tail -f /home/gene/xuchen/2.log
• 更简单的代码：

#!/bin/bash
#上一步的输出目录
export CLEAN_OUT=/home/gene/xuchen/output_gp
#设置这一步的输出目录
export WORK=/home/gene/xuchen/diffout_gp
#规范命名后的染色体模板数据
export GRCh38=/home/gene/xuchen/ratgene
#GTF数据，一般使用不带chr的
export GTF=/home/lxb/gene_template/rat/gtf/Rattus_norvegicus.Rnor_6.0.90.gtf

LIST=""
LABLES=""
for file in `ls $CLEAN_OUT`
do 
file=`echo $file  sed 's/\(.*\).sorted.bam$/\1/'`
echo $file
LIST=$LIST" "$CLEAN_OUT/$file".sorted.bam"
LABLES=$LABLES","$file
done
LABLES=`echo ${LABLES: 1}`
echo $LIST
echo $LABLES

cd $CLEAN_OUT
cuffdiff -L $LABLES -o $WORK -b $GRCh38 -p 12 -u $GTF $LIST

• 最终结果

第四步 转换成symbol的csv

gp <- read.table('~/gene/xuchen/diffout_gp/genes.fpkm_tracking', header = T)
sampleN <- stringr::str_detect(colnames(gp),'FPKM')
sampleN <- colnames(gp)[sampleN]
gp_filter <- gp[,c('gene_id', 'gene_short_name', sampleN)]
write.csv(gp_filter, file='gp.csv')

gp的内容

gp.csv的内容