昨天我們在生信技能樹解讀了 從招聘信息看一個合格的生信工程師該會哪些 ，朋友圈有一小撮人冷嘲熱諷，說辛辛苦苦學了那麼多，工資也就萬把塊錢每個月。
我就呵呵，誰讓你當年選擇了生物呢？你本來應該是去賣保險信用卡股票期貨的你，因為遇到了生信，有幸找到一份養家餬口的工作你還得寸進尺？
好不容易你才走入了社會正常運轉的一環，成為一個不可或缺的螺絲釘，對沒有背景的你來說，已經是鯉躍龍門啦！
我不知道到底多少錢的工資才符合你的期望，社會上反正多得是一夜暴富，收入顯赫的職業，你做不來，在我們生信技能樹這裡抱怨是沒有用的。但是，我仍然是可以給你秀一下，學生物信息學數據分析，也不是不可以掙錢。記住，是掙錢，合理的掙錢，不是賺錢！
話說，也許是五年前，或者是四年前，我還在寫生信菜鳥團博客的時候，分享了大量的RNA-seq分析實戰教程，非常的受歡迎。遇到一個粉絲，跟著我的教程學了好久，至少郵件問了我不下十個問題，最後還是搞不定。其實就是因為沒有學Linux，自己的Windows電腦搞虛擬機折騰起來累人。然後博士畢業需要用那個數據，確實沒辦法，就乾脆找我，開價800塊錢，我幫忙拿到表達矩陣。
公共數據首先下載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
代碼都是三五年前就寫好的：
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 這個是其中一個比較的結果，可以看出相關性很不錯，說明兩個流程的結果大致是一致的趨勢。
親愛的粉絲們，大家覺得這個項目，8個樣本的RNA-seq數據的分析，值800塊錢嗎？歡迎留言參與哦！
文末友情宣傳強烈建議你推薦給身邊的博士後以及年輕生物學PI，多一點數據認知，讓他們的科研上一個臺階：
不點讚也不打賞，為什麼呢？