我用这个技能一杯咖啡的功夫就挣了800块钱

昨天我们在生信技能树解读了 从招聘信息看一个合格的生信工程师该会哪些 ，朋友圈有一小撮人冷嘲热讽，说辛辛苦苦学了那么多，工资也就万把块钱每个月。

我就呵呵，谁让你当年选择了生物呢？你本来应该是去卖保险信用卡股票期货的你，因为遇到了生信，有幸找到一份养家糊口的工作你还得寸进尺？

好不容易你才走入了社会正常运转的一环，成为一个不可或缺的螺丝钉，对没有背景的你来说，已经是鲤跃龙门啦！

我不知道到底多少钱的工资才符合你的期望，社会上反正多得是一夜暴富，收入显赫的职业，你做不来，在我们生信技能树这里抱怨是没有用的。但是，我仍然是可以给你秀一下，学生物信息学数据分析，也不是不可以挣钱。记住，是挣钱，合理的挣钱，不是赚钱！

话说，也许是五年前，或者是四年前，我还在写生信菜鸟团博客的时候，分享了大量的RNA-seq分析实战教程，非常的受欢迎。遇到一个粉丝，跟着我的教程学了好久，至少邮件问了我不下十个问题，最后还是搞不定。其实就是因为没有学Linux，自己的Windows电脑搞虚拟机折腾起来累人。然后博士毕业需要用那个数据，确实没办法，就干脆找我，开价800块钱，我帮忙拿到表达矩阵。

公共数据

• GEO地址是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE35005

• 

• SRA地址是：https://trace.ncbi.nlm.nih.gov/Traces/sra?study=SRP010262

首先下载sra文件

内容如下：

下载sra文件

然后转为fastq文件

可以看到，都是单端测序数据，如下：

得到fastq测序数据

这些fastq文件需要根据项目信息来合并，因为一个样本测到了两个fastq文件。如下：

项目信息

接着走RNA-seq的上游分析流程

主要是拿到bam文件和表达矩阵

bam文件

流程代码，真心很简单，但是是对会的人来说，这个ngs上游分析能力，真的很简单，但是如果你并没有学过Linux的shell编程，下面代码的确是难如天书！

# for i in {0..5};do bash step1-fastp-hisat2-se.sh `pwd` config.se 6 $i;done
analysis_dir=$1
config_file=$2
number1=$3
number2=$4
conda activate rna
index=/home/yb77613/reference/index/hisat/hg38/genome
# for config.se
# sampleName and fastq files(for single end)
# SCBO-9_orgP9    /home/bladder/1.1.raw_fq/SCBO-9_orgP9.fastq.gz
cat $config_file |while read id
do
echo $id
    arr=($id)
    fq=${arr[1]}
    sample=${arr[0]}
    if((i%$number1==$number2))
    then
    ## step1: basic QC for fastq files
    if [  ! -f  ok.fastp.$sample.status ]; then
          fastp -i $fq -o 2.1.clean_fq/${sample}.fq.gz \
          --thread=4 --length_required=36 --n_base_limit=6 \
          --compression=6 -R ${sample}  \
          1>0.0.logs/log.fastp.${sample}.txt 2>&1
    fi ## end for if files
    if [ $? -eq 0 ]; then
             touch  ok.fastp.$sample.status
                else
                     echo "fastp failed" $sample
    fi
    fastqc $fq  -O 0.1.qc
    fastqc 2.1.clean_fq/${sample}.fq.gz -O 0.1.qc
    ## step2: alignment by hisat2
    if [  ! -f  ok.hisat2.$sample.status ]; then
          hisat2  -p 4  -x  $index -U  2.1.clean_fq/${sample}.fq.gz  | \
          samtools sort  -O bam  -@ 4 -o - > 3.1.hisat2_align/${sample}.se.bam
    fi ## end for if files
    if [ $? -eq 0 ]; then
             touch  ok.hisat2.$sample.status
          else
             echo "hisat2 failed" $sample
    fi

fi # check the number1 number2
    i=$((i+1))
done


对所有的bam文件，统一进行定量流程

代码都是三五年前就写好的：

gtf="/home/reference/gtf/gencode/gencode.v32.annotation.gtf"
bin_featureCounts="/home/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts";
$bin_featureCounts -T 4  -t exon -g gene_name  -a $gtf -o all.counts.id.txt  ../hisat2/*.bam 1>counts.id.log 2>&1


有了表达矩阵就走一个差异分析咯

我也不记得我分享过多少次代码，多到大家都责怪我在炒冷饭了。

RNA-seq我们在生信技能树应该是至少推出了400篇教程，而且是我们全国巡讲的标准品知识点，其中还有一个阅读量过两万的综述翻译及其细节知识点的补充：

• RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同

• RNAseq数据，下载GEO中的FPKM文件后该怎么下游分析

相信大家听完了我B站的RNA-seq分析流程后，对这个数据的应用方向都不陌生。

简化一下结果交付给对方

 

• 首先是hisat.stat.xlsx这个excel表格里面保存着本次对原始测序数据的比对情况，包括比对率，多比对情况等等。

• 然后是raw_reads_matrix.csv包含着8个样本的近5万个在GENCODE数据库记录着的基因的表达量矩阵。

• 而filter_reads_matrix.csv就是把5万个基因里面过滤掉那些低表达量基因表达量矩阵。

• 最后的rpkm.txt就是把filter_reads_matrix.csv的reads的counts转换为RPKM值，你可以跟paper做比较。

GSM860182_SG.A_correlation.pdf 这个是其中一个比较的结果，可以看出相关性很不错，说明两个流程的结果大致是一致的趋势。

• 如果你进行后续分析的话，可以直接拿rpkm.txt即可。

• 还有那个 mouse_sperm_RNA_seq_bw_files.zip里面是比对的bw文件，可以用IGV打开之后定位具体某个基因看看上面的测序深度的分布情况。

亲爱的粉丝们，大家觉得这个项目，8个样本的RNA-seq数据的分析，值800块钱吗？欢迎留言参与哦！