查看原文
其他

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

生信技能树 生信技能树 2022-06-07

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

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

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

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

话说,也许是五年前,或者是四年前,我还在写生信菜鸟团博客的时候,分享了大量的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篇教程,而且是我们全国巡讲的标准品知识点,其中还有一个阅读量过两万的综述翻译及其细节知识点的补充:

相信大家听完了我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块钱吗?欢迎留言参与哦!

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存